波动方程作为一类典型的偏微分方程, 在物理中常被用来描述物理量的时间导数与空间导数之间的关系. 在光波、声波等在介质中的传播这些实际问题的研究中, 必然受到一些外部因素(白噪声)的影响. 因此, 随机波动方程得到了广泛的研究, 具有十分重要的现实意义. 近年来, 关于随机波动方程的解的行为这一课题得到了深入研究, Walsh[1]最先将高斯白噪声驱动的随机微分方程带入人们的视野, Millet[2]等人接着证明了带非线性多项式项的随机波动方程弱解的存在性和唯一性, 对于带阻尼项的随机波动方程, Bo[3]等人证明带非线性阻尼项的随机波动方程局部解的爆破性, Jiang和Wei[4]证明了一类带记忆性的随机阻尼波动方程解的唯一存在性, Wang[5]证明了在某些特定条件下, 带阻尼项的随机波动方程的解几乎处处指数稳定. 在偏微分方程数值模拟方法中, 差分法是最常用的方法之一, 这类研究十分广泛[6-8], 包括显示差分、隐式差分、紧差分、交替方向差分等. 在本文中, 对于随机波动方程的模拟, 我们同样用差分法来实现.
在本文中我们主要考虑如下带阻尼项的随机波动方程[9]
其中$ \alpha >0 $, 表示阻尼系数, $ \Delta u\left( t, x \right) $表示$ u\left( t, x \right) $的拉普拉斯算子, $ W\left( t, x \right) $是无穷维Q-Wiener过程, $ V\left( u\left( t, x \right) \right) $为非线性项, 通常为势能函数的导数, $ \; {{u}_{0}}\left( x \right)\in {{H}_{0}}\text{=}{{L}^{2}}\left( \left[ 0, 1 \right] \right) $、$ {{v}_{0}}\left( x \right)\in {{H}_{-1}}\left( \rm{sobolev} \right) $.
令$ v\left( t, x \right)={{u}_{t}}\left( t, x \right) $, 那么(1.1)式可以写为
将$ x $区间离散化为$ \left\{ {{x}_{i}}=\frac{i}{n} \right\}_{i=0}^{n}, n\in \mathbb{N} $, 则离散形式方程为
令$ A $为一个三对角矩阵, 表示为
记$ V\left( {{u}^{\left( N \right)}} \right)={{u}^{\left( N \right)T}}A{{u}^{\left( N \right)}}+\sum\limits_{k=0}^{N}{P\left( u_{k}^{\left( N \right)} \right)} $, 其中$ P\left( u \right) $通常表示粒子的除弹性势能外的势能函数, 且这里$ {{u}^{\left( N \right)}}={{\left( u_{0}^{\left( N \right)}, u_{1}^{\left( N \right)}, \ldots , u_{N}^{\left( N \right)} \right)}^{T}} $.
则系统的哈密顿函数为
根据[10], 我们可以得到离散形式的不变分布为
上述分布在无限维空间中的极限形式为
令$ {{H}_{0}}\text{=}{{L}^{2}}\left( \left[ 0, 1 \right] \right) $为包含所有$ \left[ 0, 1 \right] $上平方可积函数的Hilbert空间, 内积为$ {{\langle \cdot \cdot \rangle }_{0}} $, 范数为$ {{\left| \cdot \right|}_{0}} $.设$ \left\{ {{e}_{i}}\left( x \right)=\frac{1}{\sqrt{2}}\sin ix, i=1, 2, \cdots \right\} $为$ {{H}_{0}} $空间的一组正交基, 特征值为$ \{ {{\lambda }_{i}}={{i}^{2}}, i=1, 2, \cdots\} $. 将$ C_{0}^{\infty }\left( \left[ 0, 1 \right] \right) $(所有具有紧支撑的无穷次可微函数)记为$ {{H}_{s}} $, 范数为$ {{\left| u \right|}_{s}}={{\left| {{\left( \rm{-}\Delta \right)}^{\frac{s}{2}}}u \right|}_{0}} $, $ s\in \mathbb{R} $.当$ s>r $, 由嵌入定理可得, $ {{H}_{s}}\subset {{H}_{r}} $.
下面对时空白噪声进行处理, 概率空间为$ \left( \Omega , \mathcal{F}, P \right) $、$ \sigma \left\{ {{F}_{t}}, t\ge 0 \right\} $满足一般性条件, 则$ {{\mathbb{R}}^{+}}\times \left[ 0, 1 \right] $上的时空白噪声$ \overset{\centerdot }{\mathop{W}}\, \left( t, x \right) $表示$ {{\Omega }^{\rm{+}}}\times {{B}_{f}}\left( {{\mathbb{R}}^{+}}\times \left[ 0, 1 \right] \right) $上的高斯随机测度$ W $的分布导数, 即$ \overset{\centerdot }{\mathop{W}}\, \left( t, x \right)\text{=}\frac{{{\partial }^{2}}W\left( t, x \right)}{\partial t\partial x} $, 且满足
根据Walsh鞅测度理论[11], 对于$ f\in {{L}^{2}}\left( \left[ 0, t \right]\times \left[ 0, 1 \right] \right) $, 我们可以通过随机积分定义一个实值连续高斯过程$ {{\left\{ {{W}_{t}}\left( f \right) \right\}}_{t\ge 0}} $, $ W\left( t, f \right)=\int_{0}^{t}{\int_{0}^{1}{f\left( s, x \right)W\left( dsdx \right)}} $.
特殊地, 当$ f=\varphi \in {{H}_{0}} $, $ W\left( t, \varphi \right)=\int_{0}^{t}{\int_{0}^{1}{\varphi \left( x \right)W\left( dsdx \right)}} $, 这是一个一维布朗运动, 方差为$ \int_{0}^{1}{{{\left| \varphi \left( x \right) \right|}^{2}}dx={{\langle \varphi , \varphi \rangle }_{0}}} $.
因此我们可以得到一个$ {{H}_{0}} $上的柱布朗运动(cylindrical Brown motion)
这里$ \beta =\left\{ {{\beta }_{i}}\left( t \right)=\int_{0}^{t}{\int_{0}^{1}{{{e}_{i}}\left( x \right)W\left( dsdx \right)}}, t>0 \right\} $是一族独立同分布的标准布朗运动. 即时空白噪声可以写为
或
因此(1.1)可以改写为
令$ H={{H}_{0}}\times {{H}_{-1}} $, 则$ H $是Hilbert空间, 内积为
记$ A=\left( \begin{matrix} 0 & I \\ \Delta & 0 \\ \end{matrix} \right), B={{\left( -\Delta \right)}^{\frac{1}{2}}} $, 根据[12]中的泛函计算, 我们得到
则随机波动方程的弱解形式为
这里$ k\left( t \right)={{B}^{-1}}\sin \left( tB \right), \overset{\centerdot }{\mathop{k}}\, \left( t \right)=\cos \left( tB \right) $.
在波动方程数值模拟方法中, 有限差分方法是最常用的方法之一, 有限差分方法是将连续微分方程的求解区域根据时空维度划分为有限个网格, 并对网格节点的值依次进行求解.
根据文献[13, 14], 首先对求解区域进行网格划分, 取两个正整数$ mn $, 规定空间步长为$ h=\frac{1}{m} $, 时间步长$ \tau =\frac{T}{n} $, $ \; {{\Omega }_{h}}=\left\{ {{x}_{i}}|{{x}_{i}}=ih, 0\le i\le m \right\} $, $ \; {{\Omega }_{t}}=\left\{ {{t}_{k}}|{{t}_{k}}=k\tau , 0\le k\le n \right\} $, 则求解区域为$ {{\Omega }_{h\tau }}={{\Omega }_{h}}\times {{\Omega }_{\tau }} $, 格点上值为$ u=\left\{ u_{i}^{k}|0\le i\le m, 0\le k\le n \right\} $.
(1) 显式差分
我们将原方程组弱化为在仅网格点上成立的随机偏微分方程, 即得到(1.1)的离散形式
将各节点处的偏导数用中心差商近似代替得到,
将以上三式代入(3.1), 得到以下差分格式
上述格式引入了虚拟越界点$ u_{i}^{-1} $, 我们假设第一式对$ k=0 $也成立, 通过代换那么$ u_{i}^{-1} $即可消去.记网比$ r={}^{\tau }/{}_{h} $, 整理(3.5)式后得到以下三层显格式:
在稳定性条件$ r\le 1 $成立时, 对不加随机项的波动方程, 这样的差分格式的解唯一存在且是收敛到精确解的[13].
(2) 隐式差分
在隐式差分法中, 我们对空间二阶偏导数取为隐格式, 即
再用二阶中心差商代替为
经过同样的网比设定以及对虚拟越界点的处理, 隐式差分方程组为
为了模拟计算的方便, 将(3.9)式写成矩阵形式
并用追赶法依次求解每个时间层上的网格信息.对于隐式差分方法而言, 无论网比r为何值, 对不加随机项的波动方程, 隐式差分格式的解唯一存在且是收敛到精确解的[13].
在数值模拟中, 我们取空间步长$ \tau \rm{=0}\rm{.01} $, 时间步长$ h=0.01 $, 总时长$ T=20 $, 波速$ a=1 $, 满足稳定性条件, 阻尼系数$ \alpha =0.5 $, $ {{u}_{0}}\left( x \right)=\sin \pi x $, $ {{u}_{1}}\left( x \right)=\sin \pi x $.
显式及隐式差分全时段波形图如下:
从上面两图中明显可以看出隐式差分法得到的波动幅度小于显式差分法, 可能是由于隐式差分法收敛速度更快导致.
选取显式及隐式差分固定间隔的10个质点, 并绘出他们的频率图如下:
可以看出上面两图靠近波段两端的质点更集中于0点附近, 而中段质点波动性更强, 符合波动方程规律.