气候与环境研究
ClimaticandEnvironmentalResearch
Vol.24No.4
Jul.2019
WangPengfei,
王鹏飞,李建平,黄刚.2019.高阶Runge-Kutta-Li算法对二维线性平流方程的计算检验[J].气候与环境研究,24(4):417-429.EnvironmentalResearch(inChinese),24(4):417-429.doi:10.3878/j.issn.1006-9585.2019.18169
LiJianping,HuangGang.2019.HighorderRunge-Kutta-Lischemetosolvetwo-dimensionallinearadvectionequations[J].Climaticand
高阶Runge-Kutta-Li算法对二维线性平流方程的
计算检验
王鹏飞1,2
李建平3,4
黄刚1
1中国科学院大气物理研究所大气科学和地球流体力学数值模拟国家重点实验室，北京1000292中国科学院大气物理研究所季风系统研究中心，北京1001903中国海洋大学物理海洋教育部重点实验室，山东青岛266003
4青岛海洋科学与技术国家实验室区域海洋动力学与数值模拟功能实验室，山东青岛266237
摘要利用高阶Li空间微分方案(Li,2005)，实现了时间积分为3～6阶Runge-Kutta-Li(RKL)格式的求解算法。
二维线性平流方程的试验结果表明：在计算稳定的条件下，各阶算法的计算误差随时间的推移基本上是线性增加的。非转动背景场的平流算例中（高斯型的初值），高阶RKL算法可以取得较好的计算效果。与3、4、5、6阶RK算法配合的Li空间差分方案有效阶数可以达到5、7、9、10阶。RK算法的阶数为5（6）阶时，总误差控制在10-7（10-8）以内。随RK阶数增加Li微分的有效阶数有增加趋势，且总误差逐渐减小。定常转速的背景场算例中（偏心的高斯型初值），当RK阶数为3时，最优空间差分阶数为10；相应的阶数为4、5、6时对应的空间最优阶为16，22，22，总计算误差可以控制在10-15～10-16。随着精度的提高，误差的绝对值减小很迅速，说明算法是非常有效的。对于圆锥型初值（定常转速的背景场），4、5、6阶RK算法和3阶算法的效果差不多。高阶算法对此类具有导数不连续点的算例，效果不如高斯初始场好，结果不能保持正定，有些地方误差出现下冲和上翘。随着空间差分精度的提高，非正定的解数量和数值减小，误差的绝对值减小，说明了算法在一定程度上是有效的，但并不适合追求极高的算法阶数。这与谱方法中的导数不连续问题有些相似，误差的产生主要源于导数的不连续性，差分类方法仅能获得与导数连续性阶数相当的算法精度。各种算例中，采用恰当的边界条件是必要的，例如旋转背景场算例，比较适合使用无穷远边界条件，否则会出现计算不稳定或无法将计算误差控制到较小的范围内。关键词Runge-Kutta-Li格式
高阶算法
二维平流方程
中图分类号
P435
文献标识码A
文章编号1006-9585(2019)04-0417-13doi:10.3878/j.issn.1006-9585.2019.18169
HighOrderRunge-Kutta-LiSchemetoSolveTwo-dimensionalLinear
AdvectionEquations
WANGPengfei1,2,LIJianping3,4,andHUANGGang1
1StateKeyLaboratoryofNumericalModelingforAtmosphericSciencesandGeophysicalFluidDynamics,InstituteofAtmosphericPhysics,ChineseAcademyofSciences,Beijing100029
收稿日期2018-12-26；网络预出版日期2019-04-21
作者简介王鹏飞，男，1973年出生，博士，高级工程师，主要从事数值模式、并行计算、非线性可预报性等方面的研究。E-mail:wpf@mail.
iap.ac.cn
资助项目国家重点研发计划2018YFA0605904，国家自然科学基金资助项目41530426、41831175、41425019，海洋局国际合作项目GASI-IPOVAI-03
FundedbyNationalKeyResearchandDevelopmentProgramofChina(Grant2018YFA0605904),NationalNaturalScienceFoundationofChina
(Grants41530426,41831175,and41425019),SOAInternationalCooperationProgramonGlobalChangeandAir-SeaInteractions(GrantGASI-IPOVAI-03)
418
气候与环境研究
ClimaticandEnvironmentalResearch
2CenterforMonsoonSystemResearch,InstituteofAtmosphericPhysics,ChineseAcademyofSciences,Beijing100190
24卷Vol.24
3KeyLaboratoryofPhysicalOceanographyInstituteforAdvancedOceanStudies,OceanUniversityofChina,Qingdao,ShandongProvince266003
4LaboratoryforRegionalOceanographyandNumericalModeling,QingdaoNationalLaboratoryforMarineScienceandTechnology,Qingdao,ShandongProvince266237
AbstractInthisstudy,aimingtotakefulladvantagesofLi’shigh-orderspatialdifferentialmethod(Li,2005),we
implementthehybridRunge-Kutta-Li(RKL)schemetosolveatwo-dimensional(2D)linearadvectionequation.Theresultsindicatethatthecomputationerrorincreasedlinearlywithtime.Theexperimentswithano-rotatebackgroundfieldoftheGaussianinitialvaluesbyRKLschemecouldobtainapreciseresult.Theeffectivespatialorderswere5,7,9,and10,correspondingtotemporalordersof3,4,5,and6,respectively.Thefifth-(sixth-)orderRunge-Kutta(RK)integrationschemewiththeninth-(tenth-)orderLi’sdifferenceschemeinspatialdirectioncontrolledtheerrorwithin10−7(10−8).TheeffectiveorderofLi’sscheme(Li,2005)tendedtoincreasewiththeincreaseintheRKorder,andthetotalerrorgraduallydecreased.
AnotherrotatedbackgroundfieldintegratedfromaneccentricGaussian-typeinitialfeaturedsimilarresults.Theeffectivespatialorderwas10whenathird-orderRKschemewasapplied,andtheyincreasedto16,22,and22whentheorderofRKschemechangedto4,5,and6,respectively.Thecomputationerrorcouldbecontrolledwithin10−15-10−16,andthepeakofGaussianinitialvalueswerewellmaintained.TheerrordecreasedsharplywhiletheorderofRKLschemeincreased,whichindicatesthattheRKLisveryeffectivetodealwithsuchproblem.
TheexperimentsofRKLschemetosolveaconeinitialcase(withrotatedbackgroundfiled)indicatethatthefourth,fifth,andsixthRKintegrationobtainedalmostthesameprecisionresultasthethird-orderRKscheme.Thehigh-orderschemewasnotaseffectiveasitwasfortheGaussianinitialconditionwhenitaddressedaproblemthathaddiscontinuousderivates.Thecomputedsolutionwasnotpositiveinthewholegrids,andinsomeplaces,theerrorwasdownward,whileitwasupwardinsomeotherplaces.Theincreaseofspatialordercouldmaketheerrorsmaller,buttheerrordescentwasnotverysharp.Thisresultsuggeststhatthehigh-orderspatialdifferenceschemehassomebenefits,butweshouldnotexpectthatanultra-high-orderschemewillleadtoanultra-highpreciseresult.Thisphenomenonrevealsthatthehigh-orderschemeislimitedbythecontinuouspropertyoftheinitialcondition,andasaresult,theerrororderisdirectlyproportionaltotheorderofderivatesoftheinitialcondition.
TheproperboundaryconditionsareimportantfortheabovecomputationcaseswhentheRKLschemeisapplied.Forinstance,inthecomputationoftherotatedbackgroundcases,thevalueoutsidethegridboxtendsto0at∞,whichisafeasibleboundarycondition.Animproperboundaryconditionmaycausethecomputationtobeunstable,ortheerrorcannotbecontrolledtoanacceptablerange.Keywords
Runge-Kutta-Lischeme,High-orderscheme,Two-dimensionaladvectionequation
1引言
数值模式和数值模拟是定量研究天气和气候变化的主要工具，但是由于数学模型中各变量之间复杂的非线性作用关系，导致模拟结果中存在着不确定性。这些不确定性的来源既包含物理参数的选取、模式的框架、观测的误差，也包括了计算误差的影响（vonNeumannandGoldstine,1947）。计算误差对数值模拟结果的影响可以从复杂的大气环流模式（王鹏飞等，2007）、耦合模式（陈显尧等,2008)运行结果看出，也可以从简单的混沌动力系统（Lietal.,2000;Liao,2008;Wangetal.,2014）、
准地转模式（Teixeiraetal.,2007）的数值试验得到验证。因此，如何采取有效的方法来控制计算误差的增长，对长时间的数值计算和精确的数值模拟至关重要。
大气和海洋流体力学模式中有许多形如∂F∂t=LF（其中L为包含变量F以及F关于空间
变量导数的算子的发展方程，t为时间），使用数值于空间和时间方向的算法精度。对这类方程计算误
模式求解这类方程组时，计算误差的大小主要取决差的研究回顾可以参考王鹏飞等（王鹏飞等，2019）的文章。一般流体问题的空间高阶算法可以参考（Tal-Ezer,1986,19;Lele,1992;季仲贞和王
4期No.4
王鹏飞等：高阶Runge-Kutta-Li算法对二维线性平流方程的计算检验
WANGPengfeietal.HighOrderRunge-Kutta-LiSchemetoSolveTwo-dimensionalLinearAdvectionEquations
419
斌,1994;MaandFu,1996;吴声昌和刘小清,1996）。大气科学领域的高阶算法可以参考季仲贞和王斌（1994）、Li（2005）、冯涛和李建平（2007）的工作。时间积分中的误差变化可以参考Lietal.（2000）、Teixeiraetal.（2007）、Wangetal.（2012）等的结果。
大气环流模式主要包括动力框架和物理过程两个部分，平流方案包含于动力框架中，并被框架所计算的背景流场所驱动，负责平流多种大气成分（也称示踪物），如水汽、云水、气溶胶等。平流过程是物理参数化和动力过程之间的重要桥梁，它的好坏直接影响到大气成分的时空分布。尽管平流方程的形式比较简单，但要正确、高效的对其进行计算仍有许多问题需要解决。大部分的平流方案建立在欧拉框架下，空间被离散为具有一定分辨率的网格，流体微团会在不同的网格之间迁移，而迁移的过程中逐渐和周围的微团混合，若梯度计算存在误差，则在这种混合的作用下会进一步减小梯度，即存在“数值扩散”。
研究表明增加算法阶数来减小时间积分时的计算误差是有效的，在冯涛和李建平（2007）研究的基础上，Wang（2017）发展了一种新的高阶Taylor-Li-Wang（TLW）算法，它的特点是可以调节时间积分方案的阶数，完成从3～30阶甚至更高阶数的算法。他们的结果显示时间积分方案的阶数大于3之后，对应的最优空间差分精度阶数可以比6阶提高很多。如算法的时间积分阶数达到5时，一维线性平流方程的算例中相匹配的空间差分阶数为10阶，误差为10-6量级。而无粘Burgers（Hopf,1950）的算例中，与5阶时间积分相匹配的空间差分阶数为16阶，误差为10-14量级。
通常，对于一些需要高精度，但又不要求极端的高精度数值解时，5、6阶的时间积分方案就够用了。相比于TLW算法，6阶以内的RK方法编程实现要简单，因此，王鹏飞等(2019)将冯涛和李建平（2007）所用的RK3算法，改进为3～6阶可调的RK方法，配合Li高阶微分公式形成Runge-Kutta-Li（简记RKL）算法。但该文仅给出了一维系统的算例，本研究将RKL算法推广到求解二维平流方程的问题中，研究RKL方法对二维平流方程的模拟效果，评估其计算性能。
2
RK方法与Li空间差分算法联合
的RKL算法
2.1
3～6阶Runge-Kutta-Li算法
RKL算法的详情可参考（王鹏飞等,2019），
这里仅作简述。RK方法是一个经典的算法，常被用于求解微分方程，许多文献都对此方法进行了介绍（Haireretal.,1993;Butcher,2008）。
以一维平流方程为例，说明RKL算法的实现，方程的形式为
∂u在求解的过程中，时间方向使用∂t+∂∂ux=0.（1）RK方法，各
阶RK方法具体如下：3阶RK方法（Heun）为
ìï
Zzn
1=,ïïZn2=z+1τf(Z1),ï
3íïï
Zzn
+2（2）3=ï3τf(Z2),
ïZn+1=zn13
î+4τf(Z1)+4τf(Z3),
4阶RK方法为
ìï
Zn
1=z,ïïïïZ=zn
+1τfï2ï
2(Z1),ïíZ=zn+1τfï
3
2(Z2),
ïïïïZï4=zn+τf(Z3),ïïî
Zn+1=zn+16τ[f(Z1)+2f(Z2)+2f(Z3)+f(Z4)],
（3）
当RK算法阶数提升时，计算步数增加，方程的个数也增加明显，为了简便书写，可将其写为系数表的形式，例如RK4写成系数表的形式为
01122
12
012，（4）
100111116336
其中，第1行表示：Z1=zn；第2行表示：Z2=
zn+1
12τf(Z1)，竖线右侧的第1列
(2
)为τf(Z1)的系420
气候与环境研究
ClimaticandEnvironmentalResearch
24卷Vol.24
数；第3行表示：Z1
3=zn+
2
τf(Z2)，竖线右侧的第1,2列(0,1
2
)分别为τf(Z1)、τf(Z2)的系数；第4
行表示：Z4=zn+τf(Z3)，竖线右侧的第1，2，3列（0，0，1）分别为τf(Z1)、τf(Z2)、τf(Z3)的系数；依此类推。第5行（横线下）表示：Zn+1=
zn+1
6
τ[f(Z1)+2f(Z2)+2f(Z3)+f(Z4)]，横线下
为计算下一步数值Zn+1，横线上的计算均为临时变
量值。
5阶RK方法（Butcher,2008）（这里仅给出系数表）为
01144111
488120012，（5）
339
4
16-3388161
-386777-1287770
321232790
90
909090
6阶RK方法（Butcher,2008）（这里仅给出系数表）为
0
1
13
32
23
0
31
111312
3-12525553515.（6）8
-444881
3111620
-24-1182101-2613343118328026013156-3919539131111413
2000
404025200
需要说明的是，3阶以上的RK方法，系数形式并不唯一，这里选用了其中一种常用的形式。
空间方向采用任意阶精度的Li（2005）算法，
计算公式为
f=1n(m)y(m)(yi)hm∑dn+1,i,jf(yj)，
（7）
j=0
此公式的精度为(n-m+1)阶，其中d(m)n+1,i,j的含义和计算方法参考Li（2005），本文选m≡1，因此，n阶精度需要（n+1）个点即可，需要注意的是，此公式中的yi位置为其在从0到n个计算点之间的相对位置，而不是u(xi)中实际网格的绝对位置。
当x方向使用N个网格点时，
∂(u)可以使用公式（7）计算，且能保证其达到n∂阶精度，使用
xn+1个格点（n+1为从N个网格点中选取出的格
点数，n≤N-1）。例如计算xi点时，如果n+1为奇数，这n+1个格点可选为(xi-n/2,⋯,xi,⋯,xi+n/2)xi-n/2对应y0，…，xi对应yn/2，…，xi+n/2对应yn），使xi位于中间点；如果n+1为偶数，这n个格点
可选为(xi-(n+1)/2,⋯,xi,⋯,xi+(n+1)/2-1)，使xi尽量靠近中间点；如果xi靠近格点的左边界或右边界，而无法将其放于中心位置时，可以选取在左侧或右侧边界上能够覆盖xi的连续n个格点即可。对于能使用周期边界条件而避免边界效应的方程，应尽量选取周期边界条件，这样可以保证xi尽量靠近计算的中间点。
当m≡1时，空间方向1阶导数的计算公式为
f1)(y1n(1)y(i)=h∑dn+1,i,jf(yj).（8）
j=0计算时，空间导数的精度为n阶精度，为了保证空间导数能达到一定的精度，至少要求使用n+1个格点。例如取N=，若n=5阶，则每次计算导数时需要6个格点，依此类推。具体的计算公式如下：
d(1)d(1)2,0,1=1，2,1,0=-1，
（9）
-j)d(1)(-1)(1a(0)n-1,i,j
n+1,i,j=
j!(n-j)!
（，当i≠j），（10）
d(1)n+1,i,i=-
j=∑n
d
(1)n+1,i,j
（，当i=j），
（11）
0,j≠i
a(0)n-1,i,j=a0(-i,⋯,k-i,⋯,n-i),(k≠i,k≠j)=(-i)⋅(⋯)⋅(k-i)⋅(⋯)⋅(n-i),(k≠i,k≠j)（12）
求解时每步逐个格点先计算空间差，后使用
（4期No.4
王鹏飞等：高阶Runge-Kutta-Li算法对二维线性平流方程的计算检验
WANGPengfeietal.HighOrderRunge-Kutta-LiSchemetoSolveTwo-dimensionalLinearAdvectionEquations
421
RK法得到下一时刻的数值，全部格点计算完毕进入下一步积分，直到目标时刻，此即RKL方法对一维线性平流方程的求解过程。2.2
二维线性平流方程的RKL算法二维线性平流方程的形式为
∂∂ut+∂∂ux+∂∂uy=0，（13）
它的求解可按2.1节中RKL算法进行，
u(1)t(xi,yj)=-∂u(xi,yj)∂u(xi,yj)≡-u(1)∂x-(1)x(xi,yj)-uy(∂xyi,yj)，
（14）
初始场为u(x,y)|t=0=e-400éx-0.5)2
+(y-0.5)2ê(úùëû
。
更一般的二维平流方程的形式（Crowley,
1968）为
∂其中，φ(x,y,∂φt)t+u∂为待∂φ求x+v∂解流∂φy=0，（15）
；u=u(x,y,t)，v=v(x,y,t)为背景流场。
对于方程（15），可改写为
φ(1)t(xi,yj)=-u(xi,yj)∂φ(∂xxi,yj,t)-v(xi,yj)∂φ(x1)(x∂i,yyj,t).
（16）
由于φ(t
i
,yj
)可以用Li方法（Li，2005）计
算出来，因此，后面的计算同一维情形，即可用RK格式求得下一时刻的数值解，且计算精度为O(hn+τM)，其中M、n分别为时间积分、空间差
分的精度阶数。
为了检验算法的性能，设计如下旋转流场试验：
单位网格1×1中，格距Δx=Δy=0.01的流场进行逆时针旋转。风场取为以ω为角速度的均匀旋转风场：u=-ω(y-yc)，v=ω(x-xc)，其中
(xc
,yc
)=(0.5,0.5)为中心位置，角速度ω=2π/T，
T=100为周期，时间步长Δt=1/100，这样经过T
时刻，旋转风场刚好回到原来的位置（经过10000个时间步长，绕中心点一周）。初始条件一般为一个定常流场。评价算法效果的辅助量为：质量守恒比：∑ϕ2
∑/ϕ/∑∑ϕ2
0，
质量分布比：ϕ0，
峰值：ϕmax，最小值：ϕmin.
初始场分布在单位方形区域中（0.5,0.65）处的二维高斯型初值：
ϕ(x,y)|-400éë(x-0.5)2
+(y-0.65)2t=0=e
êúùû
.（17）
计算时Δx=Δy=1/100，Δt=0.01，
φ(1)t=0(xi,yj)=-u0
(xxi,yj)∂φ0(φ0(x∂xi,yj)-
v0
(xi,yj)∂φ0(x∂yi,yj)，（18）
其中∂i,yj)∂φ∂x这项可用Li算法求出，0(xi,yj)∂y同样处理，这样可以得到φ1(xi,yj)，
完成一步积分，重复这个过程直到计算完成。
2.3
RKL算法的多精度实现
高阶算法在计算时还会碰到一个问题就是舍入误差，因为计算所用的时空精度阶数都可能超过10阶，例如：Wang（2017）试验中，解的数值量级为O(1)，如果仅采用双精度计算，会发现高阶算法的绝对误差有时徘徊在10-15，这恰好接
近双精度计算时的相对误差极限，所以为了观察超高阶算法的计算效果，采用多精度计算是必要
的。本研究使用了MP（Multiple-Precision）库，采用1024二进制位精度，相当于200位以上的十进制有效数字，足以区分绝对误差小到10-200的数值解。
3RKL方法的二维平流方程试验
仿照Takacs（1985）的做法，定义总误差为：E=
1(N+1)2∑N∑N(u2D-uT)，其中，uD表示数i=0j=0值解，uT表示理论解，N为平面上（x,y）方向的格点数，这里对平均误差方差取了根号，以保证总误差和原变量具有同样的量纲。3.1
试验1：RKL算法求解二维线性平流方程时计算误差随时间的变化
二维平流方程（13）和（15）属于线性的系统，可以预期，计算误差随时间的增长速度不会像非线性系统那样迅速。但是计算时间增加时，计算机进行浮点计算的次数是线性增加的，因此，误差也不会是常数值。这个试验分为两组：
422
气候与环境研究
ClimaticandEnvironmentalResearch
24卷Vol.24
第一组试验为：数值计算方程（13）（图1a、1c），x,y方向的空间步长为h=1/200，时间步长为τ=1/400，计算区域为[0,1]×[0,1]，计算时间1～10（即400～4000步），空间阶数n统一取为10。图1a蓝、红、绿、黑色分别代表时间精度为3、4、5、6阶算法的结果，其中误差值单位分别为10-3、10-5、10-6、10-7。从图1a可见，在RK算法阶数为M＝3时，计算误差在t＝1,2,3…,10时基本呈现线性增长（M＝3时，时间积分精度较低，斜率约为0.9，图1c）；对更高的空间精度阶M=4,5,6也有类似的结果，而且斜率更接近于1。
第二组试验为：数值计算方程（15）（图1b、1d），x,y方向的空间步长为h=1/200，时间步长为τ=1/100，计算区域为[0,1]×[0,1]，旋转速度为100s一周，计算时间100～1000s，间隔100s，空间阶数n统一取为20。图1b中蓝、红、绿、黑色分别代表时间精度为3、4、5、6阶算法的结果，
其中误差值单位分别为10-8、10-11、10-14、10-14。从图1b可见，在RK算法阶数为M＝3时，计算误差在t＝100,200,…,1000时基本呈现线性增长；同样对M=4,5,6也有类似的结果，而且各种时间积分精度时，斜率都非常接近于1（图1d）。
这个试验结果表明，对二维平流方程的试验，在计算稳定的前提下，计算误差随时间增加的变化趋势近似为线性的。这个性质保证了，我们若想知道方程在第5个周期或第10个周期时计算误差的大小，只需要先得到其在1个周期时的误差，然后线性倍增即可。有了这个性质，在后面的数值试验中，我们都仅对一个周期的情况进行分析，来研究误差与时间积分精度和空间差分精度的关系。3.2
试验2：RKL算法不同阶数求解二维平流系
统（定常非旋转背景场）时计算误差的变化数值求解平流方程方程（13），使用初始流函
u(x,y)|t=0=e
22
ê(x-0.5)ú-400é+(y-0.5)ùëû
数条件：
.
（19）
图1Fig.1
（a）定常非旋转背景场、（b）定常旋转的背景场计算误差随时间周期的变化；（c）定常非旋转背景场、（d）定常旋转的背景场误差(a)Thestationarynon-rotatedback-groundflowand(b)thestationaryrotatedback-groundflowerrorversusspatialdifferenceorder;the
比随时间周期的变化。蓝、红、绿、黑色分别代表时间精度为3、4、5、6阶算法的结果
abscissaisthecycleorder,theordinateistheerror;(c)and(d)aresameas(a)and(b);theordinateistheratiooferrorversuscycles.Theblue,red,green,andblackcurvesdenotethethird,fourth,fifth,andsixthtime-integrationorders,respectively
4期No.4
王鹏飞等：高阶Runge-Kutta-Li算法对二维线性平流方程的计算检验
WANGPengfeietal.HighOrderRunge-Kutta-LiSchemetoSolveTwo-dimensionalLinearAdvectionEquations
423
计算时x,y方向的空间步长为h=1/200，时间步长为τ=1/400，计算区域为[0,1]×[0,1]，计算400步，试验结果如图2，从图2b可见，在RK算法阶数M为3时，空间阶数n达到5阶时，计算误差与n大于5阶时的计算误差相差不多；M=4阶时，n达到7阶左右计算误差接*缓变化；M=5阶时，n达到9阶左右计算误差接*缓变化；M=6阶时，n达到10阶之后计算误差接*缓变
可以用来检查计算时，t时刻数值解的误差和准确度。
从图3b可见，在RK算法阶数M为3时，空间阶数n达到10阶时，计算误差与n大于10阶时的计算误差相差不多；M=4阶时，n达到16阶左右计算误差接*缓变化；M=5阶时，n达到22阶左右计算误差接*缓变化；M=6阶时，n达到22阶之后计算误差接*缓变化。需要注意的是图3b中M＝5和M＝6结果的差别不是很大，这不是由于算法性能造成的，而是因为边界条件的问题。在使用无穷远边条件时，网格点外侧的数值被设置为0，因此，在靠近分布中心点的一侧边界处会有较大的边界误差，本例中它的量级为10-22左右（e-400×0.35≈10-22，其他方向边界外也有大小不一的
2
化。这些结果说明，本试验中总计算误差受时间积分精度的影响较大，当时间方向精度大于3阶时，空间方向的格式精度可以超过6阶。3.3
试验3：RKL算法不同阶数求解二维平流系统（定常旋转背景场和Guass初值）时计算误差的变化
数值求解平流方程方程（15），使用背景流函数为
ìu=-ω(y-yc),
ïí
ïv=ωx-x,(c)î
误差），这种误差是无法通过增加算法阶数来降低的，在经历10000步积分计算后，累积到约10-16量级，这是此类旋转背景场数值求解时常见的问题
（20）（在拓展求解网格区域范围，或使分布中心点更接近于求解区域中心后，误差仍可以减小）。
表1给出了不同阶RKL方法，求解平流方程时，一个整周期后的数值结果辅助量。表中可见，n＝2时，数值解计算误差较大，在n大于等于4阶之后，误差减小得很迅速，说明了高阶空间差分的重要性。表中的数据也可以看出，各阶（M＝4,5,6表略）算法对应的最优空间差分阶数与图3b所得
其中，(xc,yc)=(0.5,0.5)为中心位置，角速度ω=2π/T，T=100为周期，时间步长τ=1/100，初始
场分布在单位方形区域中（0.50,0.65）处的二维高斯分布φ(x,y)|t=0=e
φ(x,y,t)=e
é(x-0.5)2+(y-0.65)2úù-400êëû
。
本问题在t时刻的解析解为：
é(x-0.5+0.15sinωt)2+(y-0.5-0.15cosωt)2úù-400êëû
，（21）
结果一致。
图2Fig.2
RKL方法对非旋转背景场中二维平流系统的数值试验：（a）二维Gauss型初始场；（b）计算误差随空间精度阶数的变化，横坐标为ExperimentsoftheRKLmethodfor2Dadvectionequationwithnon-rotatedback-groundflow:(a)2DGaussian-typeinitialcondition;(b)
空间精度阶数，纵坐标为误差取对数，蓝、红、绿、黑色分别代表时间精度为3、4、5、6阶
errorversusspatialdifferenceorder,wheretheabscissaisthespatialdifferenceorder,andtheordinateisthelogarithmoferror,andtheblue,red,green,andblackcurvesdenotethethird,fourth,fifth,andsixthtime-integrationorders,respectively
424
气候与环境研究
ClimaticandEnvironmentalResearch
24卷Vol.24
图3Fig.3
RKL方法对旋转背景场中二维平流系统的数值试验：（a）二维偏心Gauss型初始场；（b）计算误差随空间精度阶数的变化，横坐标ExperimentsoftheRKLmethodfor2Dadvectionequationwithrotatedback-groundflow:(a)2DeccentricGaussian-typeinitial;(b)error
为空间精度阶数，纵坐标为误差取对数，蓝、红、绿、黑色分别代表时间精度为3、4、5、6阶
versusspatialdifferenceorder,wheretheabscissaisthespatialdifferenceorder,theordinateisthelogarithmoferror,andtheblue,red,green,andblackcurvesdenotethethird,fourth,fifth,andsixthtime-integrationorders,respectively
表1试验3M＝3时2～20阶空间算法在一个周期后的数值结果
TheresultsofExpt.3forM＝3and2-20spatialordersafteronecycle
最小值φmin
-0.124604604217228-0.00016935311-0.000000761431457-0.000000019569404-0.000000000811802-0.000000000039977-0.000000000002135-0.000000000000157-0.000000000000030-0.000000000000004
质量守恒比∑φ2
Table1
阶数n2468101214161820
峰值φmax
0.8791000206055680.9921801324057090.9997254150958730.99998480446860.9999987114065180.9999998011657330.9999999167877230.9999999324514820.9999999350457870.999999935554148
0.9999999712282140.9999999674984330.99999996728730.9999999672723190.9999999672709300.9999999672707830.9999999672707610.9999999672707620.9999999672707610.999999967270759
∑φ02质量分布比∑φ
1.0000000392470051.0000000001479131.0000000000017971.0000000000000101.0000000000000060.9999999999999980.9999999999999980.9999999999999991.0000000000000000.999999999999992
∑φ0
为了进一步分析计算误差，在图4中给出了一个整周期后，解析解和数值解的图像。图4a、4c、4e分别为解析解以及n＝2、n＝4时的数值解，可见当n＝2时，可以通过观察原变量的分布，看出数值解的差异，随着计算误差的迅速减小，到n＝4时，数值解和解析解从外观上已无明显差别。此时，分析数值解和解析解的差值场就显得必要了。图4d为n＝2时误差场的分布，可见误差的量级在10-1左右，而n＝4时，通过观察误差场，可知误差的量级在10-2左右。图4b则直接给出了n＝20时误差场的情况，误差的量级在10-8左右，与图3b结果吻合，且可以观察误差的细节。
总体来看，高阶算法较好地保持了高斯初始场的最大值，但结果并不正定，有些地方误差下冲和上
翘。随着空间阶数的提高，非正定的解数量和数值减小，误差的绝对值减小，说明了算法是有效的。3.4
试验4：RKL算法不同阶数求解二维平流系
统（定常旋转背景场和圆锥初值）时计算误差的变化
数值求解平流方程方程（15），使用背景流函
ìu=-ω(y-y0),
ïí
ïv=ωx-x,(0)î
数条件：
（22）
其中(x0,y0)=(0.5,0.5)为中心位置，角速度ω=2π/T，T=100为周期，初始场分布在边长1单位的方形区
域中（0.5,0.65）处，高度为0.05，底边长为0.15的圆锥分布。时间步长τ=0.01，格距Δx=Δy=0.01，
4期No.4
王鹏飞等：高阶Runge-Kutta-Li算法对二维线性平流方程的计算检验
WANGPengfeietal.HighOrderRunge-Kutta-LiSchemetoSolveTwo-dimensionalLinearAdvectionEquations
425
图4RKL方法对旋转背景场且高斯初值的数值试验：（a）一个周期后的解析解；（b）M＝3、n＝20时，一个周期后的误差分布；（c）M＝3、n＝2时，一个周期后的数值解；（d）M＝3、n＝2时误差的分布；（e）M＝3、n＝4时，一个周期后的数值解；（f）M＝3、n＝4时误差的分布Fig.4
ExperimentsoftheRKLmethodfor2Dadvectionequationwithrotatedback-groundflowandGaussinitialcondition:(a)Theanalytical
solutionatonecycle;(b)theerrorforM＝3,n＝20atonecycle;(c)thenumericalsolutionforM＝3,n＝2atonecycle;(d)sameas(b)butforM＝3,n＝2;(e)sameas(c)butforM＝3,n＝4;(f)sameas(b)butforM＝3,n＝4
ì1
ïï30.15-φ(x,y)|t=0=íïï0,î
((x-0.5)+(y-0.65),
2
2
)0.15-0.15-
(x-0.5)+(y-0.65)>0
2
2
（23）
对于t时刻，本问题的解析解为
ì1
ï0.15-φ(x,y,t)=í3
ïî0,
(x-0.5)+(y-0.65)≤0
2
2
2
((x+0.15sinωt-0.5)+(y-0.15cosωt-0.5),
2
)φ>0φ≤0
（24）
426
表2
气候与环境研究
ClimaticandEnvironmentalResearch
试验4M＝3时2～20阶空间算法在一个周期后的数值结果
TheresultsofExpt.4forM＝3and2-20spatialordersafteronecycle
最小值φmin
-0.002800692371273-0.001137179417055-0.000680662711126-0.000576925951550-0.000527735633149-0.0004246858111-0.000497981071535-0.000420520571531-0.0004095977869-0.000401733081757
质量守恒比∑φ2
24卷Vol.24
Table2
阶数n2468101214161820
峰值φmax
0.0423187536974060.0462319880379690.0472809107352040.0477956661559950.0484805481609540.0482198506024000.0483560036281130.0483881629162220.04878799426900.0488387005057
0.9999999942096850.9999999915504280.9999999900095990.99999998234660.99999998800100.9999999874159180.9999999868552870.9999999863772550.9999999859621880.9999999855955
∑φ02质量分布比∑φ
1.0000629330945720.99996610537620.9999213297452701.0001244813312690.9997118082196211.0000544754417481.00014847829691.0001023798093741.0001302747560591.0001240684463
∑φ0
可以用来检查计算时，t时刻数值解的误差和准确度。
从图5b可见，这个圆锥型初值的误差分布曲线与图3b明显不同。首先是M＝3,4,5,6等不同阶的时间积分算法的曲线基本重合在一起，暗示此算例中时间积分精度的影响不大。其次，计算误差随空间差分阶数增加有减小的趋势，但减小的速度较慢，仅从n＝2时的10-3减小到n＝20的10-4量级。
图5c-5f可以看到，空间阶数n＝2、4、6、8时，一个整周期后数值解的分布情况。当n＝2时，可以通过观察原变量的分布，看出数值解的差异，随着n的增加，计算误差在减小，到n＝4、6、8时数值解和解析解的差别逐渐变小，这从原变量图中即可看出。
表2给出了不同空间阶数RKL方法完成试验4时，一个整周期后的数值结果统计数据。表2中可见，n＝2时，数值解计算误差较大，在n≥10阶之后，误差变化不明显，这与表1中的辅助量变化形式不同。
试验4结果说明：高阶算法对此类具有导数不连续点的算例，效果不如Guass初始场的好。数值结果不能保持正定，有些地方误差下冲和上翘。随着空间差分精度的提高，非正定的解数量和数值会减小，误差的绝对值减小，说明了算法在一定程度上是有效的，但并不适合追求极高的算法阶数。这样的结果，与谱方法中的导数不连续问题有些相似，误差的产生主要源于导数的不连续性，差分类
方法仅能获得导数连续性阶数相当的算法精度。
值得一提的是，若数值解能够具有周期边界条件，求解精度会比非周期边界条件的高，这主要是周期边界条件满足时，总可以将目标格点放置于差分公式中的中点附近，计算时误差小（如试验2）。而非周期的边界条件，在边界点以及紧靠边界点处，即使能够用n阶精度的差分格式，目标格点也会被放置于偏离n个差分点中心的位置，实际计算表明这些位置的差分精度较差。若可以近似的使用无穷远边条件（例如，边界外的值趋向于0，且很小），可以将边界外的设置为0值，以便保证相应的差分格式有足够多的有效格点来保持计算精度（如试验3和试验4）。
4结论
本文利用Li提出的高阶显式微分公式，结合Runge-Kutta时间积分方案，实现了求解含时偏微分方程的Runge-Kutta-Li高阶算法格式。对二维线性平流方程，通过比较理论解和各阶计算格式的数值解，研究了计算误差随时间积分阶数的变化情况。结果表明：在计算稳定的条件下，各算法的计算误差随时间的推移基本上是线性增加的。一个非转动背景场的平流算例中（高斯型的初值），高阶RKL算法可以取得较好的计算效果。如果时间方向的积分精度为3阶，计算误差在空间差分阶数到达一定阶数后减小的不明显。而当时间积分的阶数大于3，例如选为4、5、6阶时，此时空间差分的
4期No.4
王鹏飞等：高阶Runge-Kutta-Li算法对二维线性平流方程的计算检验
WANGPengfeietal.HighOrderRunge-Kutta-LiSchemetoSolveTwo-dimensionalLinearAdvectionEquations
427
图5RKL方法对旋转背景场且圆锥初值的数值试验：（a）二维圆锥型初始场；（b）计算误差随空间精度阶数的变化[纵坐标为误差取对
数，蓝、红、绿、黑色分别代表时间精度为3、4、5、6阶]；（c）M＝3、n＝2时，一个周期后的数值解；（d）同（c），但n＝4；（e）同（c），但n＝6；（f）同（c），但n＝8Fig.5
ExperimentsoftheRKLmethodfor2Dadvectionequationwithrotatedback-groundflowandconeinitialcondition:(a)2Dcone-typeinitial
condition;(b)errorversusspatialdifferenceorder,wheretheabscissaisthespatialdifferenceorder;theordinateisthelogarithmoferror;andtheblue,red,green,andblackcurvesdenotethethird,fourth,fifth,andsixthtime-integrationorders,respectively;(c)thenumericalsolutionsofM＝3，n＝2atonecycle;(d)sameas(c)butforn＝4;(e)sameas(c)butforn＝6;(f)sameas(c)butforn＝8
最优阶数比M＝3时的最优空间阶数要大，说明计算误差在空间差分精度达到一定阶数而出现饱和的现象是由于没有足够高的时间积分方案配合而引起的。
定常转速的背景场算例中（采用偏心的高斯型
初值），当RK算法阶数为3时，最优空间差分阶数为10；相应的M＝4,5,6时对应的空间最优阶为16、22、22，计算误差可以控制在10-15～10-16。高阶算法较好地保持了高斯初始场的最大值，随着精度的提高，误差的绝对值减小很迅速，说明算法是
428
气候与环境研究
ClimaticandEnvironmentalResearch
24卷Vol.24
非常有效的。对于圆锥型初值（定常转速的背景场），4、5、6阶RK算法和3阶算法的效果差不多。高阶算法对此类具有导数不连续点的算例，效果不如高斯型初始场的好。数值解也不能保持正定，有些地方误差下冲和上翘。随着空间差分精度的提高，非正定的解数量和数值减小，误差的绝对值减小，说明了算法在一定程度上是有效的，但并不适合追求极高的算法阶数。这与谱方法中的导数不连续问题有些相似，误差的产生主要源于导数的不连续性，差分类方法仅能获得导数连续性阶数相当的算法精度。上述算例中，采用恰当的边界条件是非常必要的，例如旋转背景场算例，比较适合使用无穷远边界条件，否则会出现计算不稳定或无法将计算误差控制到较小的范围内。
虽然可以通过减小时间步长来改进时间积分的精度，但是它的效率不如增加阶数高，这在以往的研究中已经被论证（Wangetal.,2012）。因此，相比于以往的仅增加空间差分精度的高阶算法，本文在时空两个方向提高算法的精度阶数，减少了总体计算误差。为了控制高阶算法计算过程中浮点舍入误差的影响，使用多精度计算的工具库实现了可以求解二维平流方程的RKL格式，此方法可以方便的拓展到其他类似方程的求解中。
二维平流方程还有更为复杂的背景流场，如Smolarkiewicz变形流场（Smolarkiewicz,1982;Staniforthetal.,1987）。该流场可用于研究由流函数定义的流场内标量分布的平流，反映了彼此对称的涡旋作用下平流场的运动特性。对这类较复杂的背景流场，使用欧拉框架下的算法通常无法取得特别好的计算效果，因此，能否拓展高精度算法，在拉格朗日描述（StaniforthandCôté,1991）的框架下处理这种问题，将在今后的工作中进行研究。
参考文献（References）
ButcherJC.2008.NumericalMethodsforOrdinaryDifferentialEquations(2nded.)[M].England:Wiley,463pp.
陈显尧,宋振亚,赵伟,等.2008.气候模式系统模拟结果的不确定性分析[J].海洋科学进展,26(2):119-125.
ChenXianyao,Song
Zhenya,ZhaoWei,etal.2008.Uncertainityanalysisofresultssimulatedbyclimatemodelsystem[J].AdvancesinMarineScience(inChinese),26(2):119-125.doi:10.3969/j.issn.1671-67.2008.02.001
CrowleyWP.1968.Numericaladvectionexperiments[J].Mon.Wea.Rev.,96(1):1-11.doi:10.1175/1520-0493(1968)096<0001:NAE>2.
0.CO;2
冯涛,李建平.2007.高精度迎风偏斜格式的比较与分析[J].大气科学,31(2):245-253.
FengTao,LiJianping.2007.Acomparison
andanalysisofhighorderupwind-biasedschemes[J].ChineseJournalofAtmosphericSciences(inChinese),31(2):245-253.doi:10.3878/j.issn.1006-95.2007.02.06
HairerE,NørsettSP,WannerG.1993.SolvingOrdinaryDifferentialEquationsI:NonstiffProblems[M].Berlin,Heidelberg:Springer,528pp.
HopfE.1950.Thepartialdifferentialequationut+uux＝μxx[J].Commun.PureAppl.Math.,3(3):201-230.doi:10.1002/cpa.3160030302
季仲贞,王斌.1994.一类高时间差分精度的平方守恒格式的构造及其应用检验[J].自然科学进展—国家重点实验室通讯,4(2):149-157.
JiZhongzhen,WangBin.1994.Constructionandapplicationtestofakindofhighprecisionschemewithsquare-conservation[J].ProgressinNaturalScience(inChinese),4(2):149-157.
LeleSK.1992.Compactfinitedifferenceschemeswithspectral-likeresolution[J].J.Comput.Phys.,103(1):16-42.doi:10.1016/0021-9991(92)90324-R.
LiJP.2005.Generalexplicitdifferenceformulasfornumericaldifferentiation
[J].
Journal
of
Computational
and
Applied
Mathematics,183(1):29-52.doi:10.1016/j.cam.2004.12.026LiJP,ZengQC,ChouJF.2000.Computationaluncertaintyprincipleinnonlinearordinarydifferentialequations(I)——Numericalresults[J].ScienceinChina(Ser.E),43(5):449-460.
LiaoSJ.2008.Onthereliabilityofcomputedchaoticsolutionsofnon-lineardifferentialequations[J].TellusA:DynamicMeteorologyandOceanography,61(4):550-5.doi:10.1111/j.1600-0870.2009.00402.x
MaYW,FuDX.1996.Supercompactfinitedifferencemethod(SCFDM)witharbitraryhighaccuracy[J].Comput.FluidDyn.J.,5(2):259-276.
SmolarkiewiczPK.1982.Themulti-dimensionalcrowleyadvectionscheme[J].Mon.Wea.Rev.,110(12):1968-1983.doi:10.1175/1520-0493(1982)110<1968:CAS>2.0.CO;2
StaniforthA,CôtéJ.1991.Semi-lagrangianintegrationschemesforatmosphericmodels—Areview[J].Mon.Wea.Rev.,119(9):2206-2223.doi:10.1175/1520-0493(1991)119<2206:SLISFA>2.0.CO;2StaniforthA,
Côté
J,
Pudykjewicz
J.
1987.
Comments
on
“swolarkiewicz'sdeformationalflow”[J].Mon.Wea.Rev.,115(4):4-900.doi:10.1175/1520-0493(1987)115<04:CODF>2.0.CO;2TakacsLL.1985.Atwo-stepschemefortheadvectionequationwithminimizeddissipationanddispersionerrors[J].Mon.Wea.Rev.,113(6):1050-1065.doi:10.1175/1520-0493(1985)113<1050:ATSSFT>2.0.CO;2
Tal-EzerH.1986.Spectralmethodsintimeforhyperbolicequations[J].SIAMJournalonNumericalAnalysis,23(1):11-26.doi:10.1137/0723002
Tal-EzerH.19.Spectralmethodsintimeforparabolicproblems[J].SIAMJournalonNumericalAnalysis,26(1):1-11.doi:10.1137/
4期No.4
王鹏飞等：高阶Runge-Kutta-Li算法对二维线性平流方程的计算检验
WANGPengfeietal.HighOrderRunge-Kutta-LiSchemetoSolveTwo-dimensionalLinearAdvectionEquations
429
0726001
TeixeiraJ,ReynoldsCA,JuddK.2007.Timestepsensitivityofnonlinearatmosphericmodels:Numericalconvergence,truncationerrorgrowth,andensembledesign[J].J.Atmos.Sci.,(1):175-1.doi:10.1175/JAS3824.1
vonNeumannJ,GoldstineHH.1947.Numericalinvertingofmatricesofhighorder[J].BulletinoftheAmericanMathematicalSociety,53(11):1021-1099.doi:10.1090/S0002-9904-1947-009-6
WangPF.2017.Ahigh-orderspatiotemporalprecision-matchingTaylor-Lischemefortime-dependentproblems[J].AdvancesinAtmosphericSciences,34(12):1461-1471.doi:10.1007/s00376-017-7018-1
王鹏飞,王在志,黄刚.2007.舍入误差对大气环流模式模拟结果的影响[J].大气科学,31(5):815-825.
WangPengfei,WangZaizhi,
HuangGang.2007.Theinfluenceofround-offerrorontheatmosphericgeneralcirculationmodel[J].ChineseJournalofAtmosphericSciences(inChinese),31(5):815-825.doi:10.3878/j.issn.1006-95.2007.05.06
WangPF,LiJP,LiQ.2012.Computationaluncertaintyandthe
applicationofahigh-performancemultipleprecisionschemetoobtainingthecorrectreferencesolutionofLorenzequations[J].NumericalAlgorithms,59(1):147-159.doi:10.1007/s11075-011-9481-6
WangPF,LiuY,LiJP.2014.Cleannumericalsimulationforsomechaoticsystemsusingtheparallelmultiple-precisionTaylorscheme[J].ChineseScienceBulletin,59(33):4465-4472.doi:10.1007/s11434-014-0412-5
王鹏飞,楚苹瓖,王立志,等.2019.Runge-Kutta算法与Li差分法不同阶数配合对计算精度影响研究[J].大气科学,43(1):99-106.
WangPengfei,ChuPingxiang,WangLizhi,etal.2019.A
studyontheprecisionofRunge-KuttamethodwithvariousordersofLidifferencescheme[J].ChineseJournalofAtmosphericSciences(inChinese),43(1):99-106.doi:10.3878/j.issn.1006-95.1805.17238
吴声昌,刘小清.1996.KdV方程的时间谱离散方法[J].应用数学和力学,17(4):357-362.
WuShengchang,LiuXiaoqing.1996.
SpectralmethodintimeforKDVequations[J].AppliedMathematicsandMechanics(inChinese),17(4):357-362.
因篇幅问题不能全部显示，请点此查看更多更全内容