KdV方程的一个六阶空间精度守恒差分格式

郭祯 吴明明 胡劲松

本文对含齐次边界条件的KdV方程的初边值问题进行了数值研究. 通过在时间层进行二阶精度的Crank-Nicolson差分离散、在空间层进行六阶精度的外推组合差分离散，本文建立了一个具有六阶空间精度的两层非线性差分格式.该格式能够合理地模拟原问题的两个守恒量. 然后，本文利用能量方法证明了格式的收敛性和稳定性. 数值算例验证了该方法的有效性.

KdV方程； Crank-Nicolson差分格式； 六阶精度； 守恒

O241.82A2023.031006

收稿日期： 2022-06-30

基金项目： 国家自然科学基金青年基金（11701481）；四川应用基础研究项目（2019JY0387）

作者简介： 郭祯（1997-）， 女， 硕士研究生， 主要研究方向为计算数学. E-mail： gz2080118252@163.com

通讯作者： 胡劲松. E-mail： hjs888hjs@163.com

A conservative difference scheme with 6-order spatial accuracy for the KdV equation

GUO Zhen， WU Ming-Ming， HU Jin-Song

（School of Science， Xihua University， Chengdu 610039， China）

In this paper， we propose a two-level difference scheme with 6-order spatial accuracy for the initial boundary value problem of KdV equation with homogeneous boundary condition. In this scheme， the Crank-Nicolson differential dispersion with 2-order accuracy is used in the time layer and the discretization of space layer is performed by extrapolating difference combination with 6-order accuracy. This scheme can simulate two conservative properties of the original problem reasonably. Then the convergence and stability of the scheme are proved by using the energy method. Finally， numerical examples verify the performance of the scheme.

KdV equation; Crank-Nicolson difference scheme; 6-order accuracy; Conservation

（2010 MSC 65M60）

1 引 言

KdV方程［1］

ut+αuux+uxxx=0（1）

（其中α为非零实常数）是荷兰数学家Korteweg和de Vires在研究浅水中的小振幅长波运动时提出的. 因具有无穷多个不变量，该方程在冷等离子体中的离子声波、岩浆流与管道波、泡沫液体中的声波及海洋内波等诸多研究领域有广泛应用［2］. 作为一类典型的非线性色散波动方程，KdV方程少有解析解，这就使得对其数值解法的研究具有重要价值［3-12］.

本文考虑如下KdV方程的初边值问题：

ut+αuux+uxxx=0，

x， t∈xL， xR×0， T（2）

ux， 0=u0x，x∈xL， xR（3）

uxL， t=uxR， t=0，

uxxR， t=0，t∈0， T（4）

郭 祯， 等： KdV方程的一个六阶空间精度守恒差分格式

其中u0（x）是一个已知函数. 初边值问题（2）～（4）满足如下守恒律［13］：

Qt=∫xRxLux， tdx=

∫xRxLu0x， tdx=Q0（5）

Et=u2L2=u02L2=E0（6）

其中Q（0）与E（0）均为仅与初始条件有关的常数. 该方程是奇数阶的，且没有对称性，故其数值求解有相当难度. 在已有的求解方法中，文献［9］提出了精度为Oτ2+h的两层非线性差分格式和三层线性差分格式，文献［10］构造了具有二阶精度的两层非线性差分格式.尽管这些差分格式能够模拟KdV方程的部分守恒律，但理论精度较低. 值得注意的是，为了提高数值解法的理论精度，文献［11］引入了中间函数，将方程转化为含有两个未知函数的方程组，进而提出了一个在空间层具有六阶精度的紧致差分格式. 然而，该方法需要求解两个未知函数，计算量大. 此外，文献［12］利用正交分解方法建立了一个在空间层具有六阶精度的隐式紧致差分格式，但仅验证了方法的可行性而没有对方程的守恒律进行模拟.

在本文中， 我们对初边值问题（2）～（4）建立了一个具有六阶空间精度的两层非线性差分格式，给出其收敛性和稳定性的证明. 我们的格式能够合理地模拟守恒量（5），（6）. 该格式首先在时间层进行二阶精度的Crank-Nicolson差分离散，然后在空间层进行六阶精度的外推组合差分离散. 最后我们通过数值算例验证了格式的有效性.

2 差分格式及其守恒性

记剖分区域xL， xR×0， T，时间步长为τ，tn=nτ，0≤n≤N， N=Tτ；xj=xL+jh，0≤j≤J，h=xR-xLJ为空间步长. 我们约定C是与空间步长和时间步长均无关一般常数且C>0，并定义如下记号

unj≡uxj， tn，  Unj≈uxj， tn，

Unjx=Unj+1-Unjh， Unjx-=Unj-Unj-1h，

Unjx^=Unj+1-Unj-12h， Unjx¨=Unj+2-Unj-24h，

Unjx…=Unj+3-Unj-36h， Unjt=Un+1j-Unjτ，

Un+1/2j=Un+1j+Unj2，

〈Un，Vn〉=h∑J-1j=1UnjVnj，

Un2=〈Un，Un〉， Un∞=max1≤j≤J-1Unj，

Z0h={U=Uj|U0-i=UJ+i=0， i=0， 1， 2， 3;

j=-3， -2， -1，0，1， …， J+2， J+3}.

对初边值问题（2）～（4），我们构造如下的两层非线性差分格式

Unjt+α ΨUn+1/2j+ΦUn+1/2j=0，

1≤j≤ J-1，1≤n≤N-1 （7）

U0j=u0xj，j=0， 1， …， J （8）

Un∈Z0h，n=0， 1，  …， N （9）

其中

ΨUn+1/2j=Ψ1Un+1/2j+Ψ2Un+1/2j+

Ψ3Un+1/2j，

Ψ1Un+1/2j=12{Un+1/2jUn+1/2jx^+

Un+1/2j2x^}，

Ψ2Un+1/2j=-15{Un+1/2jUn+1/2jx¨+

Un+1/2j2x¨}，

Ψ3Un+1/2j=130{Un+1/2jUn+1/2jx…+

Un+1/2j2x…}，

ΦUn+1/2j=2615Un+1/2jxx-x^-65Un+1/2jxx-x¨+

715Un+1/2jx^x^x¨.

关于格式（7）～（9）对守恒量（5）（6）的模拟效果，我们有如下结果.

定理2.1 差分格式（7）～（9）关于以下的离散能量守恒：

Qn=h∑J-1j=1Unj=Qn-1=…=Q0（10）

En=Un2=En-1=…=E0（11）

其中 n=1， 2， …， N.

证明 由条件（9）和分部求和公式［14］，我们有

h∑J-1j=1Un+1/2jUn+1/2jx^=-h∑J-1j=1Un+1/2jx^Un+1/2j.

则h∑J-1j=1Un+1/2jUn+1/2jx^=0. 同理可得

h∑J-1j=1Un+1/2jUn+1/2jx¨=0，

h∑J-1j=1Un+1/2jUn+1/2jx…=0.

又由

h∑J-1j=1Un+1/2j2x^=0，

h∑J-1j=1Un+1/2j2x¨=0，

h∑J-1j=1Un+1/2j2x…=0，

以及

h∑J-1j=1Un+1/2jxx-x^=0，  h∑J-1j=1Un+1/2jxx-x¨=0，

h∑J-1j=1Un+1/2jx^x^x¨=0，

我们有

h∑J-1j=1ΨUn+1/2j=0， h∑J-1j=1ΦUn+1/2j=0（12）

将（7）式两端乘以h后对j从1到J-1求和，注意到（12）式，我们有

h∑J-1j=1Unjt=0（13）

由Qn的定义，我们将（13）式对n递推可得（10）式.

另一方面， （7）式对2Un+1/2取内积，由（9）式和分部求和公式有

Un2t+2α〈ΨUn+1/2，Un+1/2〉+

2〈ΦUn+1/2j，Un+1/2〉=0（14）

由于

〈Ψ1Un+1/2，Un+1/2〉=

h∑J-1j=1Un+1/2j2Un+1/2jx^+

h∑J-1j=1Un+1/2j2x^Un+1/2j=

h∑J-1j=1Un+1/2j2Un+1/2jx^-

h∑J-1j=1Un+1/2j2Un+1/2jx^=0，

同理可得

〈Ψ2Un+1/2，Un+1/2〉=0，

〈Ψ3Un+1/2，Un+1/2〉=0，

即

〈ΨUn+1/2，Un+1/2〉=0（15）

又

〈Un+1/2xx-x^，Un+1/2〉=0， 〈Un+1/2xx-x¨，Un+1/2〉=0，

〈Un+1/2x^x^x¨，Un+1/2〉=0，

即

〈ΦUn+1/2，Un+1/2〉=0（16）

由En的定义，将（15）（16）式代入（14）式后对n递推可得式（11）. 证毕.

3 差分格式的收敛性和稳定性

由Taylor公式，如果函数ux足够光滑，则当h→0时我们有

ujx^=dudxx=xj+13！h2d3udx3x=xj+

15！h2d5udx5x=xj+Oh6，

ujx¨=dudxx=xj+43！h2d3udx3x=xj+ 165！h2d5udx5x=xj+Oh6，

ujx…=dudxx=xj+93！h2d3udx3x=xj+

815！h2d5udx5x=xj+Oh6，

ujxx-x^=d3udx3x=xj+14h2d5udx5x=xj+

140h4d7udx7x=xj+Oh6，

ujxx-x¨=d3udx3x=xj+34h2d5udx5x=xj+

23120h4d7udx7x=xj+Oh6，

以及

ujx^x^x¨=d3udx3x=xj+h2d5udx5x=xj+

25h4d7udx7x=xj+Oh6.

于是

32ujx^-35ujx¨+110ujx…=

dudxx=xj+Oh6（17）

2615ujxx-x^-65ujxx-x¨+715ujx^x^x¨=

d3udx3x=xj+Oh6（18）

格式（7）～（9）的截断误差定义为

rnj=unjt+Φun+1/2j+αΨun+1/2j（19）

则由Taylor公式以及（17）（18）式可知

rnj =Oτ2+h6（20）

引理3.1 U∈Z0h，恒有 Ux^2≤Ux2， Ux¨2≤Ux^2及Ux…2≤Ux^2.

证明 由于Ux2=Ux-2，由Cauchy-Schwarz不等式有

Ux^2=14h∑J-1j=1Ujx+Ujx-·

Ujx+Ujx-≤Ux2，

Ux¨2=14h∑J-1j=1Uj+1x^+Uj-1x^·

Uj+1x^+Uj-1x^≤Ux^2，

Ux…2=19h∑J-1j=1Uj+2x^+Ujx^+Uj-2x^·

Uj+2x^+Ujx^+Uj-2x^≤

h9·3·∑J-1j=1Uj+2x^2+Ujx^2+

Uj-2x^2≤Ux^2.

证毕.

定理3.2 设u（x， t）足够光滑且u0∈H2.则格式（7）～（9）的解Un以范数 · 收敛到原问题（2）～（4）的解，收敛阶为Oτ2+h6.

证明 记enj=unj-Unj. 由式（19）减去式（7）得

rnj=enjt+Φen+1/2j+αΨun+1/2j-

αΨUn+1/2j（21）

类似（16）式，有〈Φen+1/2， en+1/2〉=0. 则（21）式对2en+1/2取内积可得

〈rn， 2en+1/2〉=en2t+2α〈Ψun+1/2- ΨUn+1/2， en+1/2〉（22）

因为ux， t在有界闭区域xL， xR×0， T内具有连续三阶偏导数uxxx，所以

uL∞≤C， uxL∞≤C（23）

再由微分中值定理有

un+1/2jx=uxj+1， tn+tn+12-uxj， tn+tn+12h=

xuξj，  tn+tn+12，

其中 xj≤ξj≤xj+1. 故

un+1/2x∞≤C（24）

又由

Ψ1un+1/2j-Ψ1Un+1/2j=12un+1/2j（en+1/2j）x^+

12en+1/2j（un+1/2j）x^+（en+1/2jun+1/2j）x^-

Ψ1en+1/2j

及〈Ψ1en+1/2， en+1/2〉=0，结合（23）（24）式以及引理3.1有

〈Ψ1un+1/2-Ψ1Un+1/2， en+1/2〉=

12h∑J-1j=1［un+1/2j（en+1/2j）x^+en+1/2j（un+1/2j）x^］en+1/2j-

2h∑J-1j=1en+1/2jun+1/2j（en+1/2j）x^=

12h∑J-1j=1en+1/2j（un+1/2j）x^en+1/2j-

h∑J-1j=112h（un+1/2jen+1/2j+1en+1/2j-un+1/2jen+1/2jen+1/2j-1）=

12h∑J-1j=1en+1/2j（un+1/2j）x^en+1/2j+h2∑J-1j=1en+1/2j+1en+1/2j

（un+1/2j）x≤C（en+12+en2）（25）

同理可得

〈Ψ2un+1/2-Ψ2Un+1/2， en+1/2〉≤

C（en+12+en2）（26）

〈Ψ3un+1/2-Ψ3Un+1/2， en+1/2〉≤

C（en+12+en2）（27）

又

〈rn， 2en+1/2〉≤rn2+12（en+12+en2）（28）

将（25）～（28）式代入（22）式，整理得

en+12-en2≤τrn2+

Cτen+12+en2（29）

将（29）式由0到n-1递推求和得

en2≤e02+τ∑n-1l=0rl2+Cτ∑nl=0el2（30）

由（20）式有

τ∑n-1l=0rl2≤nτmax0≤l≤n-1rl2≤ T·O（τ2+h6）2.

再由e02=0知式等价于

en2≤O（τ2+h6）2+Cτ∑nl=0el2.

由离散的Gronwall不等式［14］，有en≤O（τ2+h6）. 证毕.

类似定理3.2的证明，我们有

定理3.3 在定理3.2的条件下，差分格式（7）～（9）的解Un以范数 · 关于初值无条件稳定.

4 数值算例

当参数α=6时，KdV方程的孤波解［15］为u（x， t）=Asech2kx-ω  t，其中A=r2，k=r2，ω=r3/22（r是常数）. 由此可知KdV方程的物理边界（渐近边界）条件满足当 x →∞时，有u（x，t）→0，t>0. 因此，我们只需取-xL0，xR0，问题（2）～（4）与KdV方程的Cauchy问题就是一致的. 在数值计算时我们取α=6，r=0.5. 问题（2）～（4）的初值函数可取为 u0x=14sech224x. 固定xL=-20，xR=40，T=32. 为验证差分格式（7）～（9）式的数值解在不同范数下的理论精度均为O（τ2+h6），定义

orderl2=log2enh，τ/e8nh2，τ8，

orderl∞=log2enh，τ∞/e8nh2，τ8∞.

对于τ和h的不同取值，数值解在不同时刻的误差及差分格式（7）～（9）精度的检验参见表1，2，差分格式（7）～（9）对不变量（5）（6）的数值模拟见表3.

从数值结果可以看出，差分格式（7）～（9）明显具有6阶空间精度，并且能够合理地模拟守恒性质（5）（6），因而是有效的.

参考文献：

［1］ Korteveg D J， de Vries G. On the change of long waves advancing in a rectangular canal and on a new type of long stationary waves ［J］. Phil Mag， 1895， 39： 422.

［2］ Crighton D G. Applications of KdV ［J］. Acta Appl Math， 1995， 39： 39.

［3］ Khattak A J. A comparative study of numerical solutions of a class of KdV equation ［J］. Appl Math Comput， 2008， 199： 425.

［4］ Benkhaldoun F， Seaid M. New finite-volume relaxation methods for the third-order differential equations ［J］. Commun Comput Phys， 2008， 4： 820.

［5］ Dutykh D， Marx C， Francesco F. Geometric numerical schemes for the KdV equation ［J］. Comp Math Math Phys+， 2013， 53： 221.

［6］ Zabusky N J， Martin D K. Interaction of "solitons" in a collisionless plasma and the recurrence of initial states ［J］. Phys Rev Lett， 1965， 15： 240.

［7］ CourtèsC， Frédéric L， Frédéric R. Error estimates of finite difference schemes for the Korteweg–de Vries equation ［J］. IMA J Numer Anal， 2020， 40： 628.

［8］ Wang H P， Wang Y S， Hu Y Y. An explicit scheme for the KdV equation ［J］. Chin Phys Lett， 2008， 25： 2335.

［9］ Shen J， Wang X， Sun Z.The conservation and convergence of two finite difference schemes for KdV equations with initial and boundary value conditions ［J］. Numer Math： Theory Me， 2020， 13： 253.

［10］ Wang X， Sun Z. A second order convergent difference scheme for the initial-boundary value problem of Korteweg-de Vires equation ［J］. Numer Meth Part D E， 2021， 37： 2873.

［11］ 赵修成， 黄浪扬. KdV方程的一个紧致差分格式 ［J］. 工程数学学报， 2015， 32： 876.

［12］ Zhang X， Zhang P. A reduced high-order compact finite difference scheme based on proper orthogonal decomposition technique for KdV equation ［J］. Appl Math Comput， 2018， 339： 535.

［13］ Miura R M， Gardner C S， Kruskal M D. Korteweg-de Vries equation and generalizations （II） ：  existence of conservation laws and constants of motion ［J］. J Math Phys， 1968， 9： 1204.

［14］ Zhou Y L. Applications of discrete functional analysis to the finite difference method ［M］. Beijing： International Academic Publishers， 1990.

［15］ Dehghan M， Shokri A. A numerical method for KdV equation using collocation and radial basis functions ［J］. Nonlinear Dynam， 2007， 50： 111.

引用本文格式：

中 文：  郭祯， 吴明明， 胡劲松. KdV方程的一个六阶空间精度守恒差分格式［J］. 四川大学学报：  自然科学版， 2023， 60：  031006.

英 文：  Guo Z， Wu M M， Hu J S. A conservative difference scheme with 6-order spatial accuracy for the KdV equation ［J］. J Sichuan Univ：  Nat Sci Ed， 2023， 60：  031006.