双曲型方程在自然科学领域有着广泛的应用背景,它反映了自然界中的振动和波动现象,比较典型的有描述弦振动的一维波动方程,以及描述声波及电磁波传播现象的方程等.我们考虑一类带小参数的二阶双曲型初边值问题

(1.1) $ \begin{equation} \left\{ \begin{array}{ll} \frac{\partial^2 u}{\partial t^2} = \varepsilon^2a(x, t)\frac{\partial^2 u}{\partial x^2}+f(x, t), \; &0\leq x \leq L, \; 0\leq t \leq T, \\ u(x, 0) = \varphi(x), \; \frac{\partial u}{\partial t}\mid_{t = 0} = \psi(x), \; &0\leq x \leq L, \\ u(0, t) = g_0(t), \; u(L, t) = g_1(t), \; &0\leq t \leq T, \end{array} \right. \end{equation} $

其中$ 0 < \varepsilon^2 \ll 1 $是小参数,待求未知量$ u(x, t) $是关于空间$ x $和时间$ t $的函数,简记为$ u $.系数$ a(x, t) $和非齐次项$ f(x, t) $满足一定的正则性条件.问题中函数$ \varphi(x), \; \psi(x), \; g_0(t), \; g_1(t) $已知.

求解双曲型方程的有限差分法方法已有很多^[1],但这些格式一般都是基于均匀网格构造的.当问题(1.1)中参数$ \varepsilon^2 $较小时,边界层上的解可能会产生剧烈的变化,此时在均匀网格上求得的数值解不能很好的逼近实际问题.而自适应移动网格方法是一种求解这类奇异摄动问题的有效方法,它通过网格的迭代变化,能使解变化较大的区域集中较多的网格点,更好地反映解的性质.

由于移动网格方法求解问题时网格是自适应变化的,关于该方法的理论研究比较困难.文献[2-13]研究了自适应移动网格方法在流体力学等问题中的应用.其中,文献[3-4]针对一阶的双曲守恒律方程提出了基于等分布原理的移动网格方法.针对一类对流扩散方程组问题,文献[5]提出了基于等弧长分布的移动网格方法,并进行了数值验证.文献[6]提出将适合问题的层适应网格作为初始网格来进行网格移动,比较了两种移动网格方法.我们将在有限差分隐格式的基础上,提出求解含小参数的双曲模型问题(1.1)的自适应移动网格方法.

本节我们构造空间移动网格上模型方程(1.1)的差分格式.设空间和时间剖分数分别为$ N $和$ M $.固定某一时刻$ t^k $,空间网格点为$ \{x_j\mid0 = x_0 < x_1 < \cdot\cdot\cdot < x_N = L\} $,在$ t^{k} $时刻记$ u_j^{k} = u(x_j, t^k) $, $ f_j^{k} = f(x_j, t^k) $, $ a_j^{k} = a(x_j, t^k) $,记时间步长为$ \tau $,空间步长$ h_j = x_j-x_{j-1} $, $ \hbar_j = (h_j+h_{j+1})/2 $, $ D^-u_j^{k} = \frac{u_{j}^{k}-u_{j-1}^{k}}{h_{j}} $, $ \delta ^2u_j^{k} = \frac{1}{\hbar_j}(D^-u_{j+1}^{k}-D^-u_j^{k}) $,我们构造移动网格上逼近问题(1.1)的隐格式如下

(2.1) $ \begin{equation} \frac{u_j^{k+1}-2u_j^{k}+u_j^{k-1}}{\tau^2} = \frac{\varepsilon^2}{2}a_j^{k}(\delta ^2u_j^{k+1}+\delta^2u_j^{k-1})+f_j^{k}, \; j = 1, 2, \cdot\cdot\cdot, N-1. \end{equation} $

该格式在均匀网格上是无条件稳定的,精度为$ O(\tau^2+h^2) $ (参见文献[1]).

令

$ \lambda1_{j}^{k+1} = \frac{\varepsilon^2\tau^2}{2\hbar_jh_{j+1}}a_j^{k+1}, \lambda2_{j}^{k+1} = \frac{\varepsilon^2\tau^2}{2\hbar_jh_{j}}a_j^{k+1}, $

$ \widetilde{\lambda1}_{j}^{k-1} = \frac{\varepsilon^2\tau^2}{2\hbar_jh_{j+1}}a_j^{k-1}, \widetilde{\lambda2}_{j}^{k-1} = \frac{\varepsilon^2\tau^2}{2\hbar_jh_{j}}a_j^{k-1}, $

$ F_j^{k} = 2u_j^{k}+\widetilde{\lambda2}_{j}^{k-1}u_{j-1}^{k-1}-(1+\widetilde{\lambda1}_{j}^{k-1}+\widetilde{\lambda2}_{j}^{k-1})u_{j}^{k-1}+\widetilde{\lambda1}_{j}^{k-1}u_{j+1}^{k-1}+\tau^2f_j^k, $

我们可将格式(2.1)化简为

(2.2) $ \begin{equation} -\lambda2_{j}^{k+1}u_{j-1}^{k+1}+(1+\lambda1_{j}^{k+1}+\lambda2_{j}^{k+1})u_j^{k+1}+(-\lambda1_{j}^{k+1})u_{j+1}^{k+1} = F_j^{k}, \; j = 1, 2, \cdot\cdot\cdot, N-1. \end{equation} $

对于初始条件和边界条件$ u(x, 0) = \varphi(x), \; u(0, t) = g_0(t), \; u(L, t) = g_1(t) $,我们用离散格式$ u_j^0 = \varphi (x_j), \; j = 0, 1, \cdots, N $, $ u_0^k = g_0(t^k), \; u_N^k = g_1(t^k), \; k = 0, 1, \cdots, M $代替.对于初始条件$ \frac{\partial u}{\partial t}\mid_{t = 0} = \psi(x) $,为了提高精度,我们引进虚网格点$ u_j^{-1} $,构造差分格式$ \frac{u_j^{1}-u_j^{-1}}{2\tau} = \psi(x_j) $.在第一层时间上针对问题(1.1)第一式采用格式

(2.3) $ \begin{equation} \frac{u_j^1-2u_j^0+u_j^{-1}}{\tau^2} = {\varepsilon^2}a_j^{0}\delta^2u_j^{0}+f_j^{0}, j = 1, 2, \cdot\cdot\cdot, N-1. \end{equation} $

消去$ u_j^{-1} $,有

(2.4) $ \begin{equation} u_j^{1} = u_j^{0}+\tau \varphi (x_j)+\frac{\varepsilon^2\tau^2}{2}a_j^{0}\delta ^2u_j^{0}+\frac{\tau^2}{2}f_j^{0}, j = 1, 2, \cdot\cdot\cdot, N-1. \end{equation} $

由$ u_j^0 = \varphi (x_j), \; j = 0, 1, \cdots, N $及(2.2)和(2.4)式,我们可以由$ u_j^k $和$ u_j^{k+1} $求出$ t^{k+2} $时刻的数值解$ u_j^{k+2}(k = 0, 1, \cdots, M-2) $.

对于模型问题(1.1),我们用空间区域的自适应移动网格法来求解.该方法通过网格迭代产生适应问题的新网格,然后通过插值等方式将旧网格上的解传递到新网格上,以保持在不同的网格分布上它们有相同的拓扑结构.我们用到的移动网格方法是基于等分布原理构建的^[3-4].等分布原理的基本思想是使解的某个误差度量在计算网格上分布得比较均匀,这种对解的误差的度量通常被称为控制函数.首先我们介绍网格生成部分的实现.

设$ x $和$ \xi $分别代表物理坐标系和计算坐标系,根据模型问题,有$ x\in [0, L] $.假定$ \xi \in [0, 1] $,我们建立$ x $和$ \xi $之间的一一坐标变换: $ x = x(\xi), \; \xi\in [0, 1] $,且$ x(0) = 0, \; x(1) = L $.由等分布原理的标准形式

(3.1) $ \begin{equation} (\omega x_\xi)_\xi = 0, \end{equation} $

我们可求得网格映射$ x = x(\xi) $.这里$ \omega $是控制函数,一般由公式

(3.2) $ \begin{equation} \omega(u(x, t)) = \sqrt{1+\rho(u_x(x, t))^2} \end{equation} $

的离散形式求得.式中$ \rho > 0 $是一个常数,在数值实验时确定.

离散网格方程(3.1),我们使用Gauss-Seidel方法来实现旧网格$ x_{j}^{[v]} $向新网格$ x_{j}^{[v+1]} $的迭代

(3.3) $ \begin{equation} w(u^{[v]}_{j+\frac{1}{2}})(x^{[v]}_{j+1}-x^{[v+1]}_j)-w(u^{[v]}_{j-\frac{1}{2}}) (x^{[v+1]}_{j}-x^{[v]}_{j-1}) = 0. \end{equation} $

从而可得新网格的计算公式

(3.4) $ \begin{equation} x^{[v+1]}_{j} = \frac{w(u^{[v]}_{j+\frac{1}{2}})x^{[v]}_{j+1}+w(u^{[v]}_{j-\frac{1}{2}})x^{[v]}_{j-1}}{w(u^{[v]}_{j+\frac{1}{2}})+w(u^{[v]}_{j-\frac{1}{2}})}. \end{equation} $

具体计算中我们通过(3.4)式不断迭代生成新网格,直至相邻两次生成的网格之间没有明显的变化时停止迭代.即给定一个误差控制量$ \epsilon $,当$ \|x^{[v+1]}-x^{[v]}\|\leq \epsilon $时将$ \{x^{[v+1]}_{j}\} $确定为最终的自适应网格.

新网格产生后,定义在旧网格上的离散解不能和新网格单元一一对应了,我们需要对新网格点上的解进行更新.我们通过以下的守恒型插值格式来获取新网格上对应的数值解.假定新网格$ \{\tilde{x}_{j}\} $和旧网格$ \{x_{j}\} $之间移动的距离很小,用$ \tilde{u}_{j+\frac{1}{2}} $和$ u_{j+\frac{1}{2}} $分别表示数值解在$ [\tilde{x}_j, \tilde{x}_{j+1}] $和$ [x_j, x_{j+1}] $的平均值,即

(3.5) $ \begin{equation} \tilde{u}_{j+\frac{1}{2}} = \frac{1}{\tilde{x}_{j+1}-\tilde{x}_j}\int _{\tilde{x}_j}^{\tilde{x}_{j+1}} u(\tilde{x}, t) {\rm d}\tilde{x}, \; \; u_{j+\frac{1}{2}} = \frac{1}{x_{j+1}-x_j}\int _{x_j}^{x_{j+1}} u(x, t) {\rm d}x. \end{equation} $

记网格移动量为$ c(x) $,即$ c(x) = x-\tilde{x} $,假定$ |c(x)|\ll 1 $,则有

(3.6) $ \begin{eqnarray} \int _{\tilde{x}_j}^{\tilde{x}_{j+1}} u(\tilde{x}, t){\rm d}\tilde{x} & = & \int _{x_j}^{x_{j+1}} u(x-c(x), t)(1-c'(x)) {\rm d}x\\ &\approx&\int _{x_j}^{x_{j+1}} (u(x, t)-c(x)u_x(x, t))(1-c'(x)) {\rm d}x\\ &\approx&\int _{x_j}^{x_{j+1}} (u(x, t)-(cu)_x) {\rm d}x\\ & = &\int _{x_j}^{x_{j+1}} u(x, t){\rm d}x-((cu)_{j+1}-(cu)_{j}), \end{eqnarray} $

这里忽略了高阶项, $ (cu)_{j} $指$ (cu) $在第$ j $个单元的值.由$ (3.5) $和$ (3.6) $式,有

(3.7) $ \begin{equation} (\tilde{x}_{j+1}-\tilde{x}_j)\tilde{u}_{j+\frac{1}{2}} = (x_{j+1}-x_j)u_{j+\frac{1}{2}}-((cu)_{j+1}-(cu)_{j}). \end{equation} $

由插值公式(3.7),新网格上的数值解可由

(3.8) $ \begin{equation} \tilde{u}_{j+\frac{1}{2}} = \frac{x_{j+1}-x_j}{\tilde{x}_{j+1}-\tilde{x}_j}{u}_{j+\frac{1}{2}} -\frac{1}{\tilde{x}_{j+1}-\tilde{x}_j}(({cu})_{j+1}-({cu})_{j}) \end{equation} $

求得,其中$ ({cu})_{j} = \frac{c_j}{2}(u_j^{+}+u_j^{-})-\frac{|c_j|}{2}(u_j^{+}-u_j^{-}) $,这里$ u_j^{+} $和$ u_j^{-} $通过构造分片线性多项式$ u_j^{\pm} = u_{j\pm \frac{1}{2}}+\frac{1}{2}(x_j-x_{j\pm 1})s_{j\pm \frac{1}{2}} $来定义.其中$ s_{j+ \frac{1}{2}} $是$ u_x $在点$ x_{j+\frac{1}{2}} $处的近似,它定义为

$ \begin{eqnarray*} s_{j+ \frac{1}{2}} = \frac{[{\rm sign}(s_{j+ \frac{1}{2}}^{+})+{\rm sign}(s_{j+ \frac{1}{2}}^{-})]|s_{j+ \frac{1}{2}}^{+}s_{j+ \frac{1}{2}}^{-}|}{|s_{j+ \frac{1}{2}}^{+}|+|s_{j+ \frac{1}{2}}^{-}|}, \; s_{j+ \frac{1}{2}}^{+} = \frac{u_{j+\frac{3}{2}}-u_{j+\frac{1}{2}}}{x_{j+\frac{3}{2}}-x_{j+\frac{1}{2}}}, \; s_{j+ \frac{1}{2}}^{-} = \frac{u_{j+\frac{1}{2}}-u_{j-\frac{1}{2}}}{x_{j+\frac{1}{2}}-x_{j-\frac{1}{2}}}. \end{eqnarray*} $

根据自适应网格的生成方法和网格迭代过程中解重新分布的方法,我们给出移动网格方法求解模型方程(1.1)的算法.具体的步骤如下.

第一步  给定求解区域$ [0, L] $的一个均匀剖分$ \{x_j\mid x_j = jL/N, j = 0, 1, \cdot\cdot\cdot, N\} $,根据已知的初始条件求出均匀网格上$ t^0 $时刻的数值解,构造控制函数并通过网格迭代公式(3.4)不断迭代产生$ t^0 $时刻的自适应网格,直到相邻两次网格之间没有明显的变化时停止迭代,得到$ t^{0} $时刻的最终网格点$ \{x_j^{0}\} $和解分布$ \{u_{j}^{0}\} $.然后由(2.4)式产生$ t^1 $时刻的解分布$ \{u_{j}^{1}\} $.此时$ k = 0 $.

第二步  在上一步产生的自适应网格上求解模型问题(1.1),利用差分格式(2.2),由已有的$ t^{k} $和$ t^{k+1} $时刻的数值解$ \{u_{j}^{k}\} $和$ \{u_{j}^{k+1}\} $求出$ t^{k+2} $时刻的数值解$ \{u_{j}^{k+2}\} $,时间向前推进一步.

第三步  用$ \{u_{j}^{k+1}\} $生成控制函数,根据网格迭代公式(3.4)将旧网格移动到新网格,然后用插值公式(3.8)得到新网格上$ t^{k+1} $和$ t^{k+2} $时刻的解分布,再利用$ t^{k+1} $时刻的解分布更新控制函数,重复移动网格和插值的过程一直到相邻两次网格之间没有明显的变化.最后得到最终网格$ \{\tilde{x}_j^{k+1}\} $上$ t^{k+1} $和$ t^{k+2} $时刻的解分布$ \{\tilde{u}_{j}^{k+1}\} $和$ \{\tilde{u}_{j}^{k+2}\} $.

第四步  如果$ t^{k+2}\leq T $,令$ k = k+1 $,转第二步;否则,算法终止.将第二步产生的数值解作为最终数值解.

下面我们根据上述算法来进行数值实验.

例4.1

$ \begin{eqnarray*} \left\{ \begin{array}{ll} \frac{\partial^2 u}{\partial t^2} = \varepsilon^2 \frac{\partial^2 u}{\partial x^2}+f(x, t), \; &0\leq x \leq 1, \; 0\leq t \leq T, \\ u(x, 0) = \frac{e^{-\frac{x}{\varepsilon}}-1}{e^{-\frac{1}{\varepsilon}}-1} +\sin(\frac{\pi }{2}x), \; &0\leq x \leq 1, \\ \frac{\partial u}{\partial t}\mid_{t = 0} = -\frac{e^{-\frac{x}{\varepsilon}}-1} {e^{-\frac{1}{\varepsilon}}-1}-\sin(\frac{\pi }{2}x), \; &0\leq x \leq 1, \\ [2mm] u(0, t) = 0, \; u(1, t) = 2e^{-t}, \; &0\leq t \leq T. \end{array} \right. \end{eqnarray*} $

这里$ f(x, t) = [-\frac{1}{e^{-\frac{1}{\varepsilon}}-1}+ (1+\frac{\varepsilon^2\pi^2}{4})\sin(\frac{\pi }{2}x)]e^{-t} $,该问题的精确解为$ u(x, t) = [\frac{e^{-\frac{x}{\varepsilon}}-1}{e^{-\frac{1}{\varepsilon}}-1}+ \sin(\frac{\pi}{2}x)]e^{-t} $.

表 1–4列出了$ \varepsilon^2 = 0.0001 $时4个固定时刻下移动网格和均匀网格上求得的数值解与精确解的$ L_2 $误差和最大模误差.表 5–6列出了$ \varepsilon^2 = 0.0025 $时4个固定时刻下移动网格上求得的数值解与精确解的$ L_2 $误差和最大模误差.实验过程借助MATLAB软件完成.表中$ N $和$ M $分别表示求解的空间网格剖分数和时间剖分数,空间方向采用自适应移动网格,时间方向采取均匀剖分. $ r $表示收敛率,由$ r = \log_2(e_{N\times M}/e_{2N\times 2M}) $计算得出,这里$ e_{N\times M} $和$ e_{2N\times 2M} $分别对应于网格剖分数为$ {N\times M} $和$ {2N\times 2M} $时的误差.算例中控制函数选取$ \omega(u(x, t)) = \sqrt{1+10000(u_x(x, t))^2} $.我们可以看到当$ \varepsilon^2 $较小时,移动网格上求得的数值解仍能保持较好的精度,明显的优于均匀网格上求得的数值解.比如表 2中, $ t = 0.5 $时,在移动网格上取网格剖分数$ {50\times 50} $时得到的最大模误差为1.952E-04.而表 4中,在均匀网格上,即使取网格剖分数$ {200\times 200} $,也无法达到这样的精度.

图 1-图 3给出了$ \varepsilon^2 = 0.0001, t = 2, N = 50 $时移动网格和均匀网格上真解及数值解的网格点分布,其中实线为真解曲线, $ \circ $为对应网格节点的数值解.图 1和图 2分别为控制函数(3.2)中$ \rho $取10和10000时的情况.显然在边界层$ x = 0 $右侧,均匀网格上获取的数值解有明显偏离真解曲线,而移动网格上获取的数值解与真解曲线有较好的吻合.并且,该方法将更多的节点集中到了解变化较大的边界层附近,而在解变化平缓的区域网格点相对稀疏,能更好的反映解的性质.从图 1和图 2可看出,控制函数中$ \rho $的选取不同,移动到边界层的网格点的数目便不同,从而最终获取的数值解的精度也会不同.因此,我们可以通过调整$ \rho $的取值来控制网格点的移动.图 4给出了$ t = 0, N = 50 $时网格的迭代过程,可以看到初始的均匀网格逐渐向边界层移动.

图 5和图 6分别给出了$ x\in[0, 1], t\in[0, 2], N\times M = 50\times 50 $时移动网格上数值解和均匀网格上真解的图像.从整体上看,数值解与真解的值是非常吻合的.两个图像的不同在于,各时刻图 5的边界层都集中了较多节点,而图 6中由于网格点的均匀分布导致边界层的图像不能很好的刻画,这也正是移动网格方法的优越性的体现.

本文针对一类二阶双曲型方程提出了自适应移动网格方法,该方法对不同时间层的空间区域采用自适应移动网格获取数值解.我们通过数值实验说明了对于解在局部存在急剧变化的二阶双曲型问题,在移动网格上求解比均匀网格上求解具有更好的逼近效果.该方法可以推广到高维的情形.