摘要:变分迭代法是一种基于变分原理, 具有高数值精度的数值格式, 目前已广泛应用于各类强非线性孤立波方程的数值求解中. 本文利用修正的变分迭代法对两类非线性方程进行研究. 该格式是对原数值方法的一种改进, 即在变分项前引入了参数h. 通过定义误差函数的离散二范数并在定义域内绘出h-曲线, 从而确定出使误差达到最小的h, 再返回原迭代过程进行求解. 同时, 参数的引入也扩大了原数值解的收敛域, 在迭代次数一定的情况下达到了数值最优. 在数值实验中, 将上述结果应用于四阶的Cahn-Hilliard方程和Benjamin-Bona-Mahoney-Burgers方程. 对于四阶的Cahn-Hilliard方程, 普通的变分迭代法绝对误差在$10^{-1}$左右, 经过修正后, 绝对误差降为$10^{-4}$, 而且修正后的方法扩大了原数值解的收敛域. 对于Benjamin-Bona-Mahony-Burgers方程, 利用带有辅助参数的变分迭代法将数值解的精度提高到$10^{-3}$, 对真解的逼近效果优于原始的变分迭代法. 此数值方法也为其他强非线性孤立波微分方程的数值求解提供了方法和参考.
关键词: 修正的变分迭代法/
非线性孤立波方程/
四阶Cahn-Hilliard方程/
Benjamin-Bona-Mahony-Burgers方程

English Abstract

全文HTML

--> --> -->

非线性科学作为现代数学物理的一个重要组成, 具有广泛的分支. 如非线性场论^[1]、玻色-爱因斯坦凝聚物^[2]、孤立波理论^[3,4]、非线性光学^[5,6]等. 而非线性发展方程本身的物理背景和孤立波解的特殊性质使得孤子理论成为当前科学发展的前沿和热点问题. 目前, 已经发展出了多种多样的用来处理各种可积方程的解析方法, 如Fokas方法^[7,8]、反散射变换^[9-11]、达布变换^[12,13]等. 最近, 光孤子的出现也引起了一大批****的关注^[14-17]. 但是由于发展型方程的强非线性, 解析方法不总是适用的. 因此数值方法的出现大大促进了孤立子理论和非线性物理的发展. 如杨建科^[18]提出的平方迭代算法来进行孤立波的稳定性分析得到了广泛的应用; Bao和Yin^[19]利用四阶时间分裂傅里叶谱方法研究了Diarc方程, 也被推广到其他非线性方程中; Antoine等^[20]使用数值方法研究了非线性薛定谔方程的动力学性质; Cockburn和Shu^[21]以及Jiang和Shu^[22]率先提出了间断有限元和加权无本质振荡(WENO)等方法, 极大丰富了数值求解偏微分方程的格式.
本文要介绍的变分迭代法与上述几种传统的数值求解方法相比, 有着求解过程灵活、收敛速度快、数值精度高等优势. He^[23]首次提出了变分迭代法, 并对变分迭代方法进行了非常清晰的讨论. 同时该方法也被推广到各类积分-微分方程^[24,25]的数值求解中. Hesameddini和Latifizadeh^[26]使用Laplace变换重构了变分迭代算法; Salkuyeh^[27]证明了变分迭代算法的收敛性. 至此, 变分迭代法的基本理论已经成形. Noor和Mohyud-Din^[28]利用He’s多项式对变分迭代法进一步改进, 该方法将变分迭代法与同伦摄动法^[29]进行结合, 大大提高了数值精度; Zayed和Rahman^[30]利用改进后的变分迭代法研究了KdV方程和Wu-Zhang方程的数值解, 说明了变分迭代法在处理高维问题上的优势.
本文第2节简要介绍原始变分迭代法的原理以及求解过程; 第3节重点介绍带有参数摄动的修正变分迭代法^[31], 该方法通过在变分项前引入参数h, 从而对原始数值算法进行改进, 提高了格式的收敛性; 第4节和第5节将上述方法应用到两类强非线性孤立波微分方程中, 其逼近效果理想, 且与问题的物理性质兼容.

考虑一般的微分方程

$ L[u(x)]+N[u(x)] = f(x), $

其中$ L[u(x)] $和$ N[u(x)] $分别代表线性项和非线性项, $ f(x) $为已知函数. 变分迭代法的主要思想是为(1)式构造一个校正泛函如下.

$ \begin{split} u_{k+1}(x) =\;& u_k(x)+\int_0^x\lambda(\eta)\big\{L[u_k(\eta)]\\ &+N[\tilde u_k(\eta)]-f(\eta)\big\}{\rm{d}}\eta, \end{split} $

其中$ \lambda(\eta) $称为广义的拉氏乘子, 可以用变分理论最佳识别. 对(2)式两边同时对$ u_k(x) $取变分δ, 即

$ \begin{split} \delta u_{k+1}(x) =\;& \delta u_k(x)+\delta\int_0^x\lambda(\eta)\big\{L[u_k(\eta)]\\ &+N[{\tilde u_k(\eta)}]-f(\eta)\big\}{\rm{d}}\eta, \end{split} $

其中${\tilde u_k(\eta)}$称为变分限制项, 即$\delta {\tilde u_k(\eta)} = 0.$
简而言之, 方程(1)的求解过程如下, 给定初值$ u_0(x) $,

$ \begin{split} u_{k+1}(x) =\;& u_k(x)+\int_0^x\lambda(\eta)\{L[u_k(\eta)]\\ &+N[\tilde u_k(\eta)]-f(\eta)\}d\eta, \end{split} $

其中$k = 0, 1, \cdots$, 进行迭代计算.
上述通过迭代算法得到解的方法便称为变分迭代法.

可以将未知的辅助参数插入到变分迭代算法(4)式中, 则算法修正如下:
给定初值$ u_0(x) $, 通过以下迭代式进行计算

$ \begin{split} u_{k+1}(x) =\;& u_k(x)+h\int_0^x\lambda(\eta)\big\{L[u_k(\eta)]\\ &+N[u_k(\eta)]-f(\eta)\big\}{\rm{d}}\eta, \end{split} $

其中$k = 0, 1, \cdots$. 通常会采取截断解$ u_k(x, t, h) $, 其中含有辅助参数h, 它保证了数值精度、收敛速度以及收敛域, 这种方法被称为带有辅助参数的变分迭代算法. 该格式的最大优点是简便, 并且能在较大的收敛域中对原方程进行近似. 通常, 可通过误差函数在给定的区域内根据二范数画出h曲线, 从而确定出h的最佳取值, 再将h 代入到迭代算法(5)式中进行求解.

考虑如下的四阶Cahn-Hilliard方程:

$ \begin{split}&u_t+\left\{b(u)\left[\gamma u_{xxx}-(\psi'(u))_x\right]\right\}_x = 0,\quad\\&~~~~~~~~~(x,t)\in (0,4)\times(0,1), \end{split}$

同时给定初值$ u(x, 0) = u_0(x). $ 其中$ \psi(u) $是关于u的光滑函数, 一般为非线性. 四阶Cahn-Hilliard方程是Cahn和Hilliard^[32]在1958年研究二元合金的分离现象时首次提出的. 随着现代科学技术的不断发展, 该方程的应用越来越广泛, 例如多相流中的界面动力学. 一些****利用经典的数值方法对四阶Cahn-Hilliard方程进行求解, 例如有限元方法^[33]、间断有限元法^[34]、多重网格法^[35]. 本文的主要工作是通过修正的变分迭代法对四阶Cahn-Hilliard方程进行数值求解, 并且与原始的格式进行了对比.
在数值实验中, 取$b(u) = 1,\; \psi(u) = \dfrac{3}{2}(1-u^2), $ $\gamma = 4, \;u_0(x) = \sin(x)-\sin(x/2).$ 该方程的精确解为$ u(x, t) = {\rm e}^{-t}\sin(x)-{\rm e}^{t/2}\sin(x/2)$. 通过原始算法, 得出下述迭代公式

$ \begin{split}\;&{u_{n + 1}}(x,t) = {u_n}(x,t) \\& - \int_0^t \left[ \frac{{\partial {u_n}(x,s)}}{{\partial s}} + 3\frac{{{\partial ^2}{u_n}(x,s)}}{{\partial {x^2}}} + 4\frac{{{\partial ^4}{u_n}(x,s)}}{{\partial {x^4}}} \right] {\rm d}s, \end{split}$

取上述初值后便可以进行计算求解. 通过迭代求解, 图1给出了$ u_{10}(x, t) $的绝对误差图像.
图 1 原始算法的误差
Figure1. Error of the original algorithm.

由图1可以看出, 在x和t较大时, 误差会偏大. 如果继续执行迭代过程, 则可以提高精度; 当n趋向无穷大时, 可以获得精确解.
利用带有辅助参数的变分迭代法, 得出以下迭代公式

$ \begin{split} & {u_{n + 1}}(x,t,h) = {u_n}(x,t,h) \\& - h\int_0^t \left[\frac{{\partial {u_n}(x,s)}}{{\partial s}} + 3\frac{{{\partial ^2}{u_n}(x,s)}}{{\partial {x^2}}} + 4\frac{{{\partial ^4}{u_n}(x,s)}}{{\partial {x^4}}} \right] {\rm d}s. \end{split} $

下面通过方程定义误差函数如下:

$ \begin{split}r_n(x,t,h) =\;& \frac{\partial u_n(x,t,h)}{\partial t}+3\frac{\partial^2 u_n(x,t,h)}{\partial x^2}\\&+4\frac{\partial^4 u_n(x,t,h)}{\partial x^4}. \end{split}$

为了找到一个合适的参数h, 图2给出了误差函数在离散网格上的二范数. 其中, 离散后的误差函数的二范数定义为
图 2 h-曲线
Figure2. h-curve

$ \frac{1}{{400}}\sum\limits_{i = 1}^{20} {\sum\limits_{j = 1}^{20} {{r_{10}}} } \left( {\frac{i}{5},\frac{j}{{20}},h} \right). $

根据图2, 选取$ h = 0.963575 $作为辅助参数, 将其代回到(8)式进行迭代求解. 结果与真解之间的绝对误差如图3所示. 可以看出, 带有辅助参数的变分迭代法不仅扩大了收敛域, 同时提高了数值精度.
图 3 带有参数的算法的误差
Figure3. Error of algorithm with parameter

图4和图5分别给出了区域内的真解与数值解图像.
图 4 精确解
Figure4. Exact solution

图 5 数值解
Figure5. Numerical solution

图6中, 固定$ x = 4 $, 给出了数值解与真解的剖面图对比.
图 6 $ x=4$时原方法和修正后方法的数值解以及真解的图像
Figure6. The numerical solution of the original method and the corrected method, as well as the image of the exact solution when $ x=4$

Benjamin-Bona-Mahony-Burgers (BBM-B)方程表示为

$\begin{split}&u_t-\gamma u_{xxt}-\alpha u_{xx}+u_{x}+\beta uu_x = 0,\quad\\&~~~~~~~~(x,t)\in (0,2)\times(0,2), \end{split}$

给定初值为$ u_0(x) $. 其中$ \alpha, \;\beta, \;\gamma $为已知常数. 上述方程由?Benjamin等^[36]提出并研究, 用于研究均匀水槽中小振幅水波的行为. 其中, $ u(x, t) $代表水平x方向上的流体速度. 作为特殊情况, 这些方程包括非线性扩散模型的Burgers方程和非线性分散介质中的长波模型的BBM方程. Karakoc和Bhowmik^[37]利用B样条的伽辽金有限元方法进行了数值模拟.
在本文的数值实验中, 取$\alpha = 1,\; \beta = 12,\;\gamma = 1, $$ u_0(x) = -\dfrac{1}{5}\dfrac{2 {\rm e}^x+{\rm e}^{2 x}}{{(1+{\rm e}^{2 x})}^2}$, 该问题有以下精确解:

$ u(x,t) = -\frac{1}{5}\frac{2{\rm e}^{x+t/5}+{\rm e}^{2(x+t/5)}}{{\left[1+{\rm e}^{2(x+t/5)}\right]}^2}. $

原始的变分格式可以利用以下迭代公式进行求解:

$ \begin{split}& u_{n+1}(x,t) =u_n(x,t)\\&~~~~ -\int_0^t \bigg[\frac{\partial u_n(x,s)}{\partial s}-\frac{\partial^2 u_n(x,s)}{\partial x^2}+\frac{\partial u_n(x,s)}{\partial x}\\&~~~~ +12u_n(x,t) +12u_n(x,t)\frac{\partial u_n(x,s)}{\partial x}\bigg] {\rm d}s.\\[-17pt]\end{split}$

图7给出了给定区域内$ u_4(x, t) $的误差图像. 在x和t偏小的情况下, 误差不会增加很快; 但是当x和t增大时, 误差会迅速增长到$ 10 $左右. 从图7不难得出以下结论: 原始的变分迭代法对于真实解的逼近效果较差, 而且在增加迭代次数的情况下, 会出现数值不稳定的情况.
图 7 原始算法的误差
Figure7. Error of the original algorithm.

故下面采用参数摄动的变分迭代法进行数值求解. 迭代公式为

$\begin{split}\;&{u_{n + 1}}(x,t,h)={u_n}(x,t,h) \\ &\; - h \int_0^t \bigg[ \frac{{\partial {u_n}(x,s)}}{{\partial s}} - \frac{{{\partial ^2}{u_n}(x,s)}}{{\partial {x^2}}} + \frac{{\partial {u_n}(x,s)}}{{\partial x}}\\& + 12{u_n}(x,t) + 12{u_n}(x,t)\frac{{\partial {u_n}(x,s)}}{{\partial x}} \bigg] {\rm d}s.\\[-15pt] \end{split}$

通过参数的引入, 本文通过误差函数定出h的最佳取值. 误差函数定义如下:

$\begin{split}& r_n(x,t,h) = \\& \frac{\partial u_n(x,s)}{\partial s}-\frac{\partial^2 u_n(x,s)}{\partial x^2} +\frac{\partial u_n(x,s)}{\partial x} \\& + 12u_n(x,t) + 12u_n(x,t)\frac{\partial u_n(x,s)}{\partial x}. \end{split}$

引入节点后, h曲线给出了h的局部最佳取值. 定义误差函数的二范数如下:

$ \frac{1}{{100}}\sum\limits_{i = 1}^{10} {\sum\limits_{j = 1}^{10} {{r_{10}}} } \left( {\frac{i}{5},\frac{j}{5},h} \right). $

从图8可以看出, 当$ h = 0.1037 $时, 整体误差达到极小值. 故在(14)式中取$ h = 0.1037 $, 进行迭代求解. 图9给出了改进得到的数值解与真解之间的误差, 不难看出带有参数摄动的数值格式大大降低了数值解的误差, 绝对误差降低到$ 10^{-3} $, 对真解的逼近达到了理想效果. 区域内不同点的真解和数值解值如表1所列. 图10和图11给出了数值解和真解的图像对比.
图 8 h-曲线
Figure8. h-curve

图 9 带有参数的算法的误差
Figure9. Error of algorithm with parameters

$ (x, t)$	原始数值解$ u_1$	修正数值解$ u_2$	精确解u
(1, 0.2)	–1.855 × 10–1	–1.855 × 10–1	–1.855 × 10–1
(1, 0.4)	–2.176 × 10–1	–1.857 × 10–1	–1.864 × 10–1
(1, 0.6)	–2.970 × 10–1	–1.860 × 10–1	–1.873 × 10–1
(1, 0.8)	–4.234 × 10–1	–1.862 × 10–1	–1.881 × 10–1
(1, 1.0)	–6.297 × 10–1	–1.863 × 10–1	–1.889 × 10–1
(1, 1.2)	–8.914 × 10–1	–1.866 × 10–1	–1.893 × 10–1
(1, 1.4)	–1.212	–1.868 × 10–1	–1.903 × 10–1
(1, 1.6)	–1.590	–1.870 × 10–1	–1.910 × 10–1
(1, 1.8)	–2.033	–1.872 × 10–1	–1.916 × 10–1
(1, 2.0)	–2.429	–1.874 × 10–1	–1.922 × 10–1

表1不同点处真解与数值解对比
Table1.Comparison of true and numerical solutions at different points

图 10 精确解
Figure10. Exact solution

图 11 数值解
Figure11. Numerical solution

本文使用具有辅助参数的变分迭代算法可以解决四阶的Cahn-Hilliard方程和BBM-B方程的数值解问题. 具有辅助参数的变分迭代算法在应用, 计算精度与效率等方面具有明显的优势. 其中辅助参数可以通过h曲线和残差函数的2范数的误差来确定. 图形和数值结果表明, 带有参数摄动的变分迭代算法适用于物理科学和工程学中出现的大部分线性和非线性问题, 优于原始的变分迭代算法.
感谢中国矿业大学数学学院杨金杰博士的讨论.