数学杂志  2022, Vol. 42 Issue (3): 259-266   PDF    
扩展功能
加入收藏夹
复制引文信息
加入引用管理器
Email Alert
RSS
本文作者相关文章
李秋赐
随机阻尼波动方程及其有限差分模拟
李秋赐    
武汉大学数学与统计学院, 湖北 武汉 430072
摘要:本文介绍了高斯时间白噪声驱动的随机阻尼波动方程, 讨论并给出了该方程的弱解形式, 并基于两种差分法对这样的解进行了数值模拟.
关键词随机波动方程    时空白噪声    差分模拟    
STOCHASTIC WAVE EQUATION WITH DAMPING AND ITS FINITE DIFFERENCE SIMULATION
LI Qiu-ci    
School of Mathematics and Statistics, Wuhan University, Wuhan 430072
Abstract: In this paper, the stochastic wave equation driven by Gaussian time white noise is introduced, and the weak solution of the equation is discussed, and the numerical simulation of the solution is carried out based on two different methods.
Keywords: stochastic wave equation     time space white noise     difference simulation    
1 引言

波动方程作为一类典型的偏微分方程, 在物理中常被用来描述物理量的时间导数与空间导数之间的关系. 在光波、声波等在介质中的传播这些实际问题的研究中, 必然受到一些外部因素(白噪声)的影响. 因此, 随机波动方程得到了广泛的研究, 具有十分重要的现实意义. 近年来, 关于随机波动方程的解的行为这一课题得到了深入研究, Walsh[1]最先将高斯白噪声驱动的随机微分方程带入人们的视野, Millet[2]等人接着证明了带非线性多项式项的随机波动方程弱解的存在性和唯一性, 对于带阻尼项的随机波动方程, Bo[3]等人证明带非线性阻尼项的随机波动方程局部解的爆破性, Jiang和Wei[4]证明了一类带记忆性的随机阻尼波动方程解的唯一存在性, Wang[5]证明了在某些特定条件下, 带阻尼项的随机波动方程的解几乎处处指数稳定. 在偏微分方程数值模拟方法中, 差分法是最常用的方法之一, 这类研究十分广泛[6-8], 包括显示差分、隐式差分、紧差分、交替方向差分等. 在本文中, 对于随机波动方程的模拟, 我们同样用差分法来实现.

在本文中我们主要考虑如下带阻尼项的随机波动方程[9]

$ \begin{equation} \left\{ \begin{aligned} \frac{{{\partial }^{2}}u\left( t, x \right)}{\partial {{t}^{2}}}&=\Delta u\left( t, x \right)-\alpha \frac{\partial u\left( t, x \right)}{\partial t}+V\left( u\left( t, x \right) \right)+\overset{\centerdot }{\mathop{W}}\, \left( t, x \right), x\in \left( 0, 1 \right), 0\le t\le T \\ u\left( 0, x \right)&={{u}_{0}}\left( x \right), \frac{\partial u}{\partial t}\left( 0, x \right)={{v}_{0}}\left( x \right), x\in \left[ 0, 1 \right] \\ u\left( t, 0 \right)&=u\left( t, 1 \right)=0 \end{aligned} \right. \end{equation} $ (1.1)

其中$ \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)式可以写为

$ \begin{equation} \left\{ \begin{aligned} du\left( t, x \right)&=v\left( t, x \right)dt \\ dv\left( t, x \right)&=\Delta u\left( t, x \right)dt+V\left( u\left( t, x \right) \right)dt-\alpha v\left( t, x \right)dt+\overset{\centerdot }{\mathop{W}}\, \left( t, x \right)dt, 0\le t\le T \\ u\left( 0, x \right)&={{u}_{0}}\left( x \right), v\left( 0, x \right)={{v}_{0}}\left( x \right), x\in \left[ 0, 1 \right] \\ u\left( t, 0 \right)&=u\left( t, 1 \right)=0. \end{aligned} \right. \end{equation} $ (1.2)

$ x $区间离散化为$ \left\{ {{x}_{i}}=\frac{i}{n} \right\}_{i=0}^{n}, n\in \mathbb{N} $, 则离散形式方程为

$ \begin{equation} \begin{aligned} du\left( t, \frac{k}{n} \right)=&v\left( t, \frac{k}{n} \right)dt \\ dv\left( t, \frac{k}{n} \right)=&\left\{ {{n}^{2}}\left[ u\left( t, \frac{k+1}{n} \right)+u\left( t, \frac{k-1}{n} \right)-2u\left( t, \frac{k}{n} \right) \right]+V\left( u\left( t, \frac{k}{n} \right) \right)-\alpha v\left( t, \frac{k}{n} \right) \right\} \\ &dt+\sqrt{n}dW\left( t, \frac{k}{n} \right). \end{aligned} \end{equation} $ (1.3)

$ A $为一个三对角矩阵, 表示为

$ \begin{equation} A=\left( \begin{matrix} -1 & 1 & 0 & \cdots & 0 \\ 1 & -2 & 1 & 0 & \cdots \\ \vdots & {} & {} & {} & {} \\ 0 & \cdots & 1 & -2 & 1 \\ 0 & \cdots & 0 & 1 & -1 \\ \end{matrix} \right). \end{equation} $ (1.4)

$ 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}} $.

则系统的哈密顿函数为

$ \begin{equation} H\left( u, v \right)=\frac{1}{2}{{\left| v \right|}^{2}}+V\left( u \right). \end{equation} $ (1.5)

根据[10], 我们可以得到离散形式的不变分布为

$ \begin{equation} {{\pi }_{N}}\left( dudv \right)\text{=}\frac{1}{{{Z}_{N}}}\exp \left\{ -\frac{2\alpha }{N}H\left( u, v \right) \right\}\prod\limits_{k=0}^{N}{d{{u}_{k}}d{{v}_{k}}}. \end{equation} $ (1.6)

上述分布在无限维空间中的极限形式为

$ \begin{equation} \pi \left( dudv \right)=\frac{1}{Z}\exp \left\{ -2c\int_{0}^{1}{{{\left| v\left( x \right) \right|}^{2}}dx} \right\}dv\times \exp \left\{ -2c\int_{0}^{1}{\left[ {{\left| \overset{\centerdot }{\mathop{u}}\, \left( x \right) \right|}^{2}}+P\left( u\left( x \right) \right) \right]dx} \right\}du. \end{equation} $ (1.7)
2 随机波动方程的弱解

$ {{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} $, 且满足

$ \begin{equation} E\left( W\left( t, x \right)W\left( s, y \right) \right)=sx, 0\le s\le t<\infty , 0\le x\le y<1. \end{equation} $ (2.1)

根据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)

$ \begin{equation} W\left( t \right)=\sum\limits_{i=1}^{\infty }{{{\beta }_{i}}\left( t \right){{e}_{i}}}. \end{equation} $ (2.2)

这里$ \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\} $是一族独立同分布的标准布朗运动. 即时空白噪声可以写为

$ \begin{equation} \overset{\centerdot }{\mathop{W}}\, \left( t, x \right)=\sum\limits_{i=1}^{\infty }{{{\overset{\centerdot }{\mathop{\beta }}\, }_{i}}\left( t \right){{e}_{i}}\left( x \right)}, \text{for}\left( t, x \right)\in {{\mathbb{R}}^{+}}\times \left[ 0, 1 \right], \end{equation} $ (2.3)

$ \begin{equation} \overset{\centerdot }{\mathop{W}}\, \left( t, x \right)=\sum\limits_{i=1}^{\infty }{{{\overset{\centerdot }{\mathop{\beta }}\, }_{i}}\left( t \right){{e}_{i}}}. \end{equation} $ (2.4)

因此(1.1)可以改写为

$ \begin{equation} \left\{ \begin{aligned} \frac{{{\partial }^{2}}{{u}^{\varepsilon }}\left( t \right)}{\partial {{t}^{2}}}& =\Delta {{u}^{\varepsilon }}\left( t \right)-\alpha \frac{\partial {{u}^{\varepsilon }}\left( t \right)}{\partial t}+V\left( u\left( t, x \right) \right)+\sqrt{\varepsilon }\overset{\centerdot }{\mathop{W}}\, \left( t \right), 0\le t\le T \\ {{u}^{\varepsilon }}\left( 0 \right)& ={{u}_{0}}, \frac{\partial {{u}^{\varepsilon }}}{\partial t}\left( 0 \right)={{v}_{0}} \\ {{u}^{\varepsilon }}\left( t, 0 \right)& ={{u}^{\varepsilon }}\left( t, 1 \right)=0. \end{aligned} \right. \end{equation} $ (2.5)

$ \begin{equation} \left\{ \begin{aligned} {{u}^{\varepsilon }}\left( t \right)&={{u}_{0}}+\int_{0}^{t}{{{v}^{\varepsilon }}\left( s \right)ds} \\ {{v}^{\varepsilon }}\left( t \right)&={{v}_{0}}+\int_{0}^{t}{\left( \Delta {{u}^{\varepsilon }}\left( s \right)+V\left( s, {{u}^{\varepsilon }}\left( s \right) \right) \right)ds+\sqrt{\varepsilon }\int_{0}^{t}{dW\left( t \right)}}. \end{aligned} \right. \end{equation} $ (2.6)

$ H={{H}_{0}}\times {{H}_{-1}} $, 则$ H $是Hilbert空间, 内积为

$ \begin{equation} \langle \left( \begin{aligned} & {{u}_{1}} \\ & {{v}_{1}} \\ \end{aligned} \right), \left( \begin{aligned} & {{u}_{2}} \\ & {{v}_{2}} \\ \end{aligned} \right) \rangle ={{\langle {{u}_{1}}, {{u}_{2}} \rangle }_{0}}+{{\langle {{v}_{1}}, {{v}_{2}} \rangle }_{-1}}. \end{equation} $ (2.7)

$ A=\left( \begin{matrix} 0 & I \\ \Delta & 0 \\ \end{matrix} \right), B={{\left( -\Delta \right)}^{\frac{1}{2}}} $, 根据[12]中的泛函计算, 我们得到

$ \begin{equation} {{e}^{At}}=\left( \begin{matrix} \cos \left( tB \right) & {{B}^{-1}}\sin \left( tB \right) \\ -B\sin \left( tB \right) & \cos \left( tB \right) \\ \end{matrix} \right). \end{equation} $ (2.8)

则随机波动方程的弱解形式为

$ \begin{equation} {{u}^{\varepsilon }}\left( t \right)=\overset{\centerdot }{\mathop{k}}\, \left( t \right){{u}_{0}}+k\left( t \right){{v}_{0}}+\int_{0}^{t}{k\left( t-s \right)V\left( s, {{u}^{\varepsilon }}\left( s \right) \right)ds}+\sqrt{\varepsilon }\int_{0}^{t}{k\left( t-s \right)dW\left( s \right)}, \end{equation} $ (2.9)

这里$ k\left( t \right)={{B}^{-1}}\sin \left( tB \right), \overset{\centerdot }{\mathop{k}}\, \left( t \right)=\cos \left( tB \right) $.

3 差分格式的建立

在波动方程数值模拟方法中, 有限差分方法是最常用的方法之一, 有限差分方法是将连续微分方程的求解区域根据时空维度划分为有限个网格, 并对网格节点的值依次进行求解.

根据文献[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)的离散形式

$ \begin{equation} \left\{ \begin{aligned} & {{\left. \frac{{{\partial }^{2}}u}{\partial {{t}^{2}}} \right|}_{\left( t{}_{k}, {{x}_{i}} \right)}}-{{\left. \frac{{{\partial }^{2}}u}{\partial {{x}^{2}}} \right|}_{\left( t{}_{k}, {{x}_{i}} \right)}}=-\alpha {{\left. \frac{\partial u}{\partial t} \right|}_{\left( {{t}_{k}}, {{x}_{i}} \right)}}-V\left( u\left( t, x \right) \right)+\overset{\centerdot }{\mathop{W}}\, \left( t, x \right), 0<i<m, 0<k\le n \\ & u\left( 0, {{x}_{i}} \right)={{u}_{0}}\left( {{x}_{i}} \right), {{\left. \frac{\partial u}{\partial t} \right|}_{\left( 0, {{x}_{i}} \right)}}={{u}_{1}}\left( {{x}_{i}} \right), 0\le i\le m \\ & u\left( {{x}_{0}}, {{t}_{k}} \right)=u\left( {{x}_{m}}, {{t}_{k}} \right)=0 . \end{aligned} \right. \end{equation} $ (3.1)

将各节点处的偏导数用中心差商近似代替得到,

$ \begin{equation} {{\left. \frac{{{\partial }^{2}}u}{\partial {{t}^{2}}} \right|}_{\left( t{}_{k}, {{x}_{i}} \right)}}\approx \frac{u\left( {{t}_{k+1}}, {{x}_{i}} \right)-2u\left( {{t}_{k}}, {{x}_{i}} \right)+u\left( {{t}_{k-1}}, {{x}_{i}} \right)}{{{\tau }^{2}}} \end{equation} $ (3.2)
$ \begin{equation} {{\left. \frac{{{\partial }^{2}}u}{\partial {{x}^{2}}} \right|}_{\left( t{}_{k}, {{x}_{i}} \right)}}\approx \frac{u\left( {{t}_{k}}, {{x}_{i+1}} \right)-2u\left( {{t}_{k}}, {{x}_{i}} \right)+u\left( {{t}_{k}}, {{x}_{i-1}} \right)}{{{h}^{2}}} \end{equation} $ (3.3)
$ \begin{equation} {{\left. \frac{\partial u}{\partial t} \right|}_{\left( {{t}_{k}}, {{x}_{i}} \right)}}\approx \frac{u\left( {{t}_{k+1}}, {{x}_{i}} \right)-u\left( {{t}_{k-1}}, {{x}_{i}} \right)}{2\tau }. \end{equation} $ (3.4)

将以上三式代入(3.1), 得到以下差分格式

$ \begin{equation} \left\{ \begin{aligned} & \frac{u_{i}^{k+1}-2u_{i}^{k}+u_{i}^{k-1}}{{{\tau }^{2}}}-\frac{u_{i+1}^{k}-2u_{i}^{k}+u_{i-1}^{k}}{{{h}^{2}}}=-\alpha \frac{u_{i}^{k+1}-u_{i}^{k-1}}{2\tau }-V\left( u_{i}^{k} \right)+\overset{\centerdot }{\mathop{W}}\, \left( {{t}_{k}}, {{x}_{i}} \right)\\ &0<i<m, 0<k\le n\rm{-}1 \\ & u_{i}^{0}={{u}_{0}}\left( {{x}_{i}} \right), \frac{u_{i}^{1}-u_{i}^{-1}}{2\tau }={{u}_{1}}\left( {{x}_{i}} \right), 0\le i\le m \\ & u\left( {{x}_{0}}, {{t}_{k}} \right)=u\left( {{x}_{m}}, {{t}_{k}} \right)=0. \end{aligned} \right. \end{equation} $ (3.5)

上述格式引入了虚拟越界点$ u_{i}^{-1} $, 我们假设第一式对$ k=0 $也成立, 通过代换那么$ u_{i}^{-1} $即可消去.记网比$ r={}^{\tau }/{}_{h} $, 整理(3.5)式后得到以下三层显格式:

$ \begin{equation} \left\{ \begin{aligned} u_{i}^{k+1}=&\left( \frac{1}{1\rm{+}{}^{\alpha \tau }/{}_{2}} \right)\left( {{r}^{2}}\left( u_{i-1}^{k}+u_{i+1}^{k} \right) +2\left( 1-{{r}^{2}} \right)u_{i}^{k}-\left( 1-\frac{\alpha \tau }{2} \right)u_{i}^{k-1}-{{\tau }^{2}}V\left( u_{i}^{k} \right){{\tau }^{2}}\right.\\ &\left.+\overset{\centerdot }{\mathop{W}}\, \left( {{t}_{k}}, {{x}_{i}} \right) \right), \qquad 1\le i\le m-1, 1\le k\le n-1, \\ u_{i}^{0}=&{{u}_{0}}\left( {{x}_{i}} \right), \qquad 0\le i\le m \\ u_{i}^{1}=&\frac{1}{2}\left( {{r}^{2}}\left( u_{i-1}^{0}+u_{i+1}^{0} \right) +2\tau \left( 1-\frac{\alpha \tau }{2} \right){{u}_{1}}\left( {{x}_{i}} \right) +2\left( 1-{{r}^{2}} \right)u_{i}^{0}-{{\tau }^{2}}V\left( u_{i}^{0} \right)\right.\\ &\left. +{{\tau }^{2}}{{\partial }_{t}}W\left( 0, {{x}_{i}} \right) \right), \qquad 1\le i\le m-1 \\ u_{0}^{k}=&u_{m}^{k}=0, \qquad 1\le k\le n . \end{aligned} \right. \end{equation} $ (3.6)

在稳定性条件$ r\le 1 $成立时, 对不加随机项的波动方程, 这样的差分格式的解唯一存在且是收敛到精确解的[13].

(2) 隐式差分

在隐式差分法中, 我们对空间二阶偏导数取为隐格式, 即

$ \begin{equation} {{\left. \frac{{{\partial }^{2}}u}{\partial {{x}^{2}}} \right|}_{\left( {{t}_{k}}, {{x}_{i}} \right)}}=\frac{1}{2}\left( {{\left. \frac{{{\partial }^{2}}u}{\partial {{x}^{2}}} \right|}_{\left( {{t}_{k-1}}, {{x}_{i}} \right)}}+{{\left. \frac{{{\partial }^{2}}u}{\partial {{x}^{2}}} \right|}_{\left( {{t}_{k+1}}, {{x}_{i}} \right)}} \right). \end{equation} $ (3.7)

再用二阶中心差商代替为

$ \begin{equation} \begin{aligned} {{\left. \frac{{{\partial }^{2}}u}{\partial {{x}^{2}}} \right|}_{\left( {{t}_{k}}, {{x}_{i}} \right)}} \approx \frac{1}{2}\left( \frac{u\left( {{t}_{k-1}}, {{x}_{i+1}} \right)-2u\left( {{t}_{k-1}}, {{x}_{i}} \right)+u\left( {{t}_{k-1}}, {{x}_{i-1}} \right)}{{{h}^{2}}}\right.\\ \left. +\frac{u\left( {{t}_{k+1}}, {{x}_{i+1}} \right)-2u\left( {{t}_{k+1}}, {{x}_{i}} \right)+u\left( {{t}_{k+1}}, {{x}_{i-1}} \right)}{{{h}^{2}}} \right). \end{aligned} \end{equation} $ (3.8)

经过同样的网比设定以及对虚拟越界点的处理, 隐式差分方程组为

$ \begin{equation} \left\{ \begin{aligned} & -\frac{{{r}^{2}}}{2}u_{i-1}^{k+1}+\left( 1+{{r}^{2}}+\frac{\alpha \tau }{2} \right)u_{i}^{k+1}-\frac{{{r}^{2}}}{2}u_{i+1}^{k+1}\\ &=2u_{i}^{k}+\frac{{{r}^{2}}}{2}\left( u_{i-1}^{k-1}+u_{i+1}^{k-1} \right)-\left( 1+{{r}^{2}}-\frac{\alpha \tau }{2} \right)u_{i}^{k-1}-{{\tau }^{2}}V\left( u_{i}^{k} \right)+{{\tau }^{2}}\overset{\centerdot }{\mathop{W}}\, \left( {{t}_{k}}, {{x}_{i}} \right) \\ &\qquad 1\le i\le m-1, 1\le k\le n-1, \\ & u_{i}^{0}={{u}_{0}}\left( {{x}_{i}} \right), 0\le i\le m, \\ & u_{i}^{1}=\frac{1}{2}\left( {{r}^{2}}\left( u_{i-1}^{0}+u_{i+1}^{0} \right)+2\tau \left( 1-\frac{\alpha \tau }{2} \right){{u}_{1}}\left( {{x}_{i}} \right)+2\left( 1-{{r}^{2}} \right)u_{i}^{0}-{{\tau }^{2}}V\left( u_{i}^{0} \right)+{{\tau }^{2}}\overset{\centerdot }{\mathop{W}}\, \left( 0, {{x}_{i}} \right) \right), \\ &\qquad 1\le i\le m-1, \\ & u_{0}^{k}=u_{m}^{k}=0, 1\le k\le n. \\ \end{aligned} \right. \end{equation} $ (3.9)

为了模拟计算的方便, 将(3.9)式写成矩阵形式

$ \begin{equation} \begin{aligned} & \left( \begin{matrix} 1+{{r}^{2}}+\frac{\alpha \tau }{2} & -\frac{{{r}^{2}}}{2} & {} & 0 & {} \\ -\frac{{{r}^{2}}}{2} & 1+{{r}^{2}}+\frac{\alpha \tau }{2} & -\frac{{{r}^{2}}}{2} & {} & {} \\ {} & {} & \ddots & {} & {} \\ {} & {} & -\frac{{{r}^{2}}}{2} & 1+{{r}^{2}}+\frac{\alpha \tau }{2} & -\frac{{{r}^{2}}}{2} \\ {} & 0 & {} & {} & 1+{{r}^{2}}+\frac{\alpha \tau }{2} \\ \end{matrix} \right)\left( \begin{matrix} u_{1}^{k+1} \\ u_{2}^{k+1} \\ \vdots \\ u_{m-1}^{k+1} \\ u_{m-1}^{k+1} \\ \end{matrix} \right)= \\ & \left( \begin{matrix} 2u_{1}^{k}+\frac{{{r}^{2}}}{2}\left( u_{0}^{k-1}+u_{2}^{k-1} \right)-\left( 1+{{r}^{2}}-\frac{\alpha \tau }{2} \right)u_{1}^{k-1}-{{\tau }^{2}}V\left( u_{1}^{k} \right)+{{\tau }^{2}}\overset{\centerdot }{\mathop{W}}\, \left( {{t}_{k}}, {{x}_{1}} \right)+\frac{{{r}^{2}}}{2}u_{0}^{k+1} \\ 2u_{2}^{k}+\frac{{{r}^{2}}}{2}\left( u_{1}^{k-1}+u_{3}^{k-1} \right)-\left( 1+{{r}^{2}}-\frac{\alpha \tau }{2} \right)u_{2}^{k-1}-{{\tau }^{2}}V\left( u_{2}^{k} \right)+{{\tau }^{2}}\overset{\centerdot }{\mathop{W}}\, \left( {{t}_{k}}, {{x}_{2}} \right) \\ \vdots \\ 2u_{m-2}^{k}+\frac{{{r}^{2}}}{2}\left( u_{m-3}^{k-1}+u_{m-1}^{k-1} \right)-\left( 1+{{r}^{2}}-\frac{\alpha \tau }{2} \right)u_{m-2}^{k-1}-{{\tau }^{2}}V\left( u_{m-2}^{k} \right)+{{\tau }^{2}}\overset{\centerdot }{\mathop{W}}\, \left( {{t}_{k}}, {{x}_{m-2}} \right) \\ 2u_{m-1}^{k}+\frac{{{r}^{2}}}{2}\left( u_{m-2}^{k-1}+u_{m}^{k-1} \right)-\left( 1+{{r}^{2}}-\frac{\alpha \tau }{2} \right)u_{m-1}^{k-1}-{{\tau }^{2}}V\left( u_{m-1}^{k} \right)+{{\tau }^{2}}\overset{\centerdot }{\mathop{W}}\, \left( {{t}_{k}}, {{x}_{m-1}} \right)+\frac{{{r}^{2}}}{2}u_{m}^{k+1} \\ \end{matrix} \right) \\ \end{aligned} \end{equation} $ (3.10)

并用追赶法依次求解每个时间层上的网格信息.对于隐式差分方法而言, 无论网比r为何值, 对不加随机项的波动方程, 隐式差分格式的解唯一存在且是收敛到精确解的[13].

4 数值模拟节

在数值模拟中, 我们取空间步长$ \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 $.

显式及隐式差分全时段波形图如下:

图 1 显式差分全时段波形图

图 2 隐式差分全时段波形图

从上面两图中明显可以看出隐式差分法得到的波动幅度小于显式差分法, 可能是由于隐式差分法收敛速度更快导致.

选取显式及隐式差分固定间隔的10个质点, 并绘出他们的频率图如下:

图 3 显式差分间隔固定的十个质点的频率图 (第一排从左到右依次为10到50, 第二排从左到右依次为60到100)

图 4 隐式差分间隔固定的十个质点的频率图 (第一排从左到右依次为10到50, 第二排从左到右依次为60到100)

可以看出上面两图靠近波段两端的质点更集中于0点附近, 而中段质点波动性更强, 符合波动方程规律.

参考文献
[1] Keda N, Watanabe S. Stochastic differential equations and diffusion processes[M]. Elsevier, 2014.
[2] Annie M, Pierre-Luc M. On a nonlinear stochastic wave equation in the plane: existence and uniqueness of the solution[J]. The Annals of Applied Probability, 2001, 11(3): 922–951.
[3] Bo Lijun, Tang Dan, Wang Yongjin. Explosive solutions of stochastic wave equations with damping on R d[J]. Journal of Differential Equations, 2007, 244(1): 170–187.
[4] Wei Tingting, Jiang Yiming. Stochastic wave equations with memory[J]. Chinese Annals of Mathematics(Series B), 2010, 31(03): 329–342. DOI:10.1007/s11401-009-0170-x
[5] 王苏鑫. 带阻尼的随机波动方程的解的长时间行为(英文)[J]. 南开大学学报(自然科学版), 2019, 52(3): 74–82.
[6] 孙成禹, 李胜军, 倪长宽, 张玉华. 波动方程变网格步长有限差分数值模拟[J]. 石油物探, 2008, 47(2): 123–140. DOI:10.3969/j.issn.1000-1441.2008.02.003
[7] 刘建康, 武贝贝. 一维边界阻尼波动方程指数稳定的半离散有限差分一致逼近格式[J]. 应用数学学报, 2018, 41(06): 832–845.
[8] 李晓旺. 两类阻尼波动方程的有限差分格式[D]. 山西大学, 2016.
[9] 张岩责编; 梁飞, 刘杰. 非线性随机波动方程[M]. 徐州: 中国矿业大学出版社. 2020.
[10] Wu Liming. Large and moderate deviations and exponential convergence for stochastic damping Hamiltonian systems[J]. Stochastic Processes and their Applications, 2001, 91(2): 204–238.
[11] John B W. An introduction to stochastic partial differential equations[M]. Springer, 1986.
[12] Michal R and Barry S. Methods of Modern mathematical physics. 2. Fourier analysis, self-adjointness, volume2[M]. Elsevier, 1975.
[13] 华冬英, 李祥贵. 微分方程的数值解法与程序实现[M]. 北京: 电子工业出版社, 2016.
[14] 余德浩, 汤华中. 微分方程数值解法(第2版.)[M]. 北京: 科学出版社, 2018.