龙格库塔格式的计算精度推导

引子

为活跃氛围，Maple技术交互研讨群（QQ836204107）于20年春节期间特地举办一场符号推导比赛，内容如下：

初描述

本文数值计算中提到的增量函数，本质上不出泰勒展式之右。在精度阶数推导中，主要思路仍是泰勒展开，只是展开项中带有也精度限制。纯粹暴力展开会计算过多不必要项，例如以下展开式中只需最高次数为5的项：
$${\left(a_n{\Delta t}^5+a_n{\Delta t}^4+a_n{\Delta t}^3+a_n{\Delta t}^2+a_n\Delta t+a_n\right)}^5$$
然而在Maple中只使用expand+coeff提取指定项系数，将面对全部展开的巨大计算负担。同时，仅人为手动计算规避冗余将遭遇繁琐计算。因此，使用Maple之前有必要进行理论上的计算优化。在证明完毕之后，感兴趣朋友进一步研究新格式的系数设置。

正文

（一）格式介绍

常微分方程数值求解中，除龙格库塔之外，存在不少同级不同阶或者不同级不同阶的计算格式，不同之处在于增量函数中的不同系数，相同之处在于将增量函数泰勒展开后都指向确定精度。以下展示3-5阶格式对应的系数a,b,c，部分涉及到自适应步长带有第四中系数cs.

restart:
#1 is Bogacki_Shampine3(2) methed
wholechoice[1]:=[0,1/2,3/4,1],[0,[1/2],[0,3/4],[2/9,1/3,4/9]],[2/9,1/3,4/9,0],[7/24,1/4,1/3,1/8]:

#2 is classic Runge_Kutta4 methed
wholechoice[2]:=[0,1/2,1/2,1],[0,[1/2],[0,1/2],[0,0,1]],[1/6,2/6,2/6,1/6],[]:

#3 is Kutta4 methed
wholechoice[3]:=[0,1/3,2/3,1],[0,[1/3],[-1/3,1],[1,-1,1]],[1/8,3/8,3/8,1/8],[]:

#4 is Fehlberg4(5) methed
wholechoice[4]:=[0,1/4,3/8,12/13,1,1/2],[0,[1/4],[3/32,9/32],[1932/2197,-7200/2197,7296/2197],[439/216,-8,3680/513,-845/4104],[-8/27,2,-3544/2565,1859/4104,-11/40]],[25/216,0,1408/2565,2197/4104,-1/5,0],[16/135,0,6656/12825,28561/56430,-9/50,2/55]:

#5 is Dormand_Prince5(4) methed
wholechoice[5]:=[0,1/5,3/10,4/5,8/9,1,1],[0,[1/5],[3/40,9/40],[44/45,-56/15,32/9],[19372/6561,-25360/2187,64448/6561,-212/729],[9017/3168,-355/33,46732/5247,49/176,-5103/18656],[35/384,0,500/1113,125/192,-2187/6784,11/84]],[35/384,0,500/1113,125/192,-2187/6784,11/84,0],[5179/57600,0,7571/16695,393/640,-92097/339200,187/2100,1/40]:

（二）计算优化

在核心部分中，增量函数中最重要的是每个k，如果证明计算精度为p，我们只需将它在(t,u)处展开到0~p次项。进而可以控制所有k皆展到0~p次项。以下是计算的优化理论部分与相关代码，测试时以7级5阶的Dormand_Prince格式为例。
$$\begin{gathered} k=f(t+a\Delta t,u+B\Delta t),B=\sum _{j=1}^p{b}_j(\Delta t{)}^{j-1} \\ k=f(t,u)+\sum _{n=1}^p\frac{1}{n!}{(a\Delta t\frac{\partial }{\partial t}+B\Delta t\frac{\partial }{\partial u})}^n(f)(t,u)+o{(\Delta t)}^p \\ =f(t,u)+\sum _{n=1}^p\frac{{(\Delta t)}^n}{n!}\sum _{k=0}^n{C}_n^k{a}^k{B}^{n-k}\frac{{\partial }^{n}f(t,u)}{\partial t^k\partial u^{n-k}}+o{(\Delta t)}^p \\ =\left(f\left(t,u\right)+\sum _{n=1}^p\left(\frac{{\left(\Delta t\right)}^n}{n!}\sum _{k=0}^n\left({C_n}^k{a}^k\left(\frac{{\partial }^n}{\partial u^{n-k}\partial t^k}f\left(t,u\right)\right){\left(\sum _{j=1}^{p-n+1}{b}_j{\left(\Delta t\right)}^{j-1}\right)}^{n-k}\right)\right)+o{(\Delta t)}^p\right) \end{gathered}$$

#阶数与格式选择
p,methedchoice:=5,5:
#计算优化部分的Maple实现，特别要注意学会unapply这个函数的使用
expand_series:=f(t,u)+add(dt^n/n!*add(binomial(n,k)*a^k*D[1$k,2$(n-k)](f)(t,u)*(add(b[j]*dt^(j-1),j=1..p-n+1))^(n-k),k=0..n),n=1..p):
coeffs_series:=coeff~(expand_series,dt,[seq](0..p)):
Strong_Taylor_series_expand:=unapply(coeffs_series,[a,b]):
#获取格式下系数
a,b,c,cs:=wholechoice[methedchoice]:
#计算格式级数
stage:=numelems(wholechoice[methedchoice][1]):
#k1赋值
k[1]:=[f(t,u),0$p]:
#计算剩余k
for i from 2 to stage do
	k[i]:=Strong_Taylor_series_expand(a[i],add(b[i][ha]*k[ha],ha=1..i-1));
end do:
#计算增量函数
test_coeff:=add(c[ha]*k[ha],ha=1..stage):

（三）一元泰勒展式

在“核心部分”中，为得到关于u的2~p阶表达式，我们不断对初始方程不断求导。高阶导数需要前几阶导数计算值，因此需要不断代换更新值。

ode[1]:=diff(u(t),t)=f(t,u(t)):
taylor_coeff[1]:=f(t,u):
subscombination:=ode[1]:
for i from 2 to p+1 do
    ode[i]:=diff(u(t),t$i)=subs([subscombination],rhs(diff(ode[i-1],t)));
    taylor_coeff[i]:=subs(u(t)=u,rhs(ode[i]))/i!;
    #print(%);
    subscombination:=subscombination,ode[i];
end do:
standard_coeff:=[seq](taylor_coeff[i],i=1..p+1):

（四）最终测试

不出意外，对于p阶格式的精度证明，standard_coeff与test_coeff做差结果是前p项均为0。在结果表示中，D为微分算子,D_{n个1，m个2}表示对第一个变量求n阶导且对第二个变量求m阶导。注意第p项表达式极为复杂，完全查看请在PC端浏览。

simplify(standard_coeff-test_coeff)

$[0,0,0,0,0,\frac{\left(D_{2,2,2,2,2}\right)\left(f\right)\left(t,u\right){\left(f\left(t,u\right)\right)}^5}{648000}+\frac{\left(10\left(D_{2,2,2}\right)\left(f\right)\left(t,u\right)\left(D_{2,2}\right)\left(f\right)\left(t,u\right)+5\left(D_{2,2,2,2}\right)\left(f\right)\left(t,u\right)D_2\left(f\right)\left(t,u\right)+5\left(D_{1,2,2,2,2}\right)\left(f\right)\left(t,u\right)\right){\left(f\left(t,u\right)\right)}^4}{648000}+\frac{\left(-15\left(D_{2,2,2}\right)\left(f\right)\left(t,u\right){\left(D_2\left(f\right)\left(t,u\right)\right)}^2+\left(-10{\left(\left(D_{2,2}\right)\left(f\right)\left(t,u\right)\right)}^2+10\left(D_{1,2,2,2}\right)\left(f\right)\left(t,u\right)\right)D_2\left(f\right)\left(t,u\right)+20\left(D_{1,2,2}\right)\left(f\right)\left(t,u\right)\left(D_{2,2}\right)\left(f\right)\left(t,u\right)+20\left(D_{2,2,2}\right)\left(f\right)\left(t,u\right)\left(D_{1,2}\right)\left(f\right)\left(t,u\right)+10\left(D_{2,2,2,2}\right)\left(f\right)\left(t,u\right)D_1\left(f\right)\left(t,u\right)+10\left(D_{1,1,2,2,2}\right)\left(f\right)\left(t,u\right)\right){\left(f\left(t,u\right)\right)}^3}{648000}+\frac{\left(165{\left(D_2\left(f\right)\left(t,u\right)\right)}^3\left(D_{2,2}\right)\left(f\right)\left(t,u\right)-45\left(D_{1,2,2}\right)\left(f\right)\left(t,u\right){\left(D_2\left(f\right)\left(t,u\right)\right)}^2-40D_2\left(f\right)\left(t,u\right)\left(D_{1,2}\right)\left(f\right)\left(t,u\right)\left(D_{2,2}\right)\left(f\right)\left(t,u\right)+\left(10{\left(\left(D_{2,2}\right)\left(f\right)\left(t,u\right)\right)}^2+30\left(D_{1,2,2,2}\right)\left(f\right)\left(t,u\right)\right)D_1\left(f\right)\left(t,u\right)+40\left(D_{1,2,2}\right)\left(f\right)\left(t,u\right)\left(D_{1,2}\right)\left(f\right)\left(t,u\right)+10\left(D_{2,2,2}\right)\left(f\right)\left(t,u\right)\left(D_{1,1}\right)\left(f\right)\left(t,u\right)+10\left(D_{2,2}\right)\left(f\right)\left(t,u\right)\left(D_{1,1,2}\right)\left(f\right)\left(t,u\right)+10\left(D_{1,1,1,2,2}\right)\left(f\right)\left(t,u\right)\right){\left(f\left(t,u\right)\right)}^2}{648000}+\frac{\left(-180{\left(D_2\left(f\right)\left(t,u\right)\right)}^5+180{\left(D_2\left(f\right)\left(t,u\right)\right)}^3\left(D_{1,2}\right)\left(f\right)\left(t,u\right)+\left(150D_1\left(f\right)\left(t,u\right)\left(D_{2,2}\right)\left(f\right)\left(t,u\right)-30\left(D_{1,1,2}\right)\left(f\right)\left(t,u\right)\right){\left(D_2\left(f\right)\left(t,u\right)\right)}^2+\left(-30\left(D_{1,2,2}\right)\left(f\right)\left(t,u\right)D_1\left(f\right)\left(t,u\right)-10\left(D_{1,1}\right)\left(f\right)\left(t,u\right)\left(D_{2,2}\right)\left(f\right)\left(t,u\right)-40{\left(\left(D_{1,2}\right)\left(f\right)\left(t,u\right)\right)}^2-10\left(D_{1,1,1,2}\right)\left(f\right)\left(t,u\right)\right)D_2\left(f\right)\left(t,u\right)+15\left(D_{2,2,2}\right)\left(f\right)\left(t,u\right){\left(D_1\left(f\right)\left(t,u\right)\right)}^2+\left(20\left(D_{1,2}\right)\left(f\right)\left(t,u\right)\left(D_{2,2}\right)\left(f\right)\left(t,u\right)+30\left(D_{1,1,2,2}\right)\left(f\right)\left(t,u\right)\right)D_1\left(f\right)\left(t,u\right)+20\left(D_{1,2,2}\right)\left(f\right)\left(t,u\right)\left(D_{1,1}\right)\left(f\right)\left(t,u\right)+20\left(D_{1,2}\right)\left(f\right)\left(t,u\right)\left(D_{1,1,2}\right)\left(f\right)\left(t,u\right)+5\left(D_{1,1,1,1,2}\right)\left(f\right)\left(t,u\right)\right)f\left(t,u\right)}{648000}-\frac{{\left(D_2\left(f\right)\left(t,u\right)\right)}^4{D}_1\left(f\right)\left(t,u\right)}{3600}+\frac{{\left(D_2\left(f\right)\left(t,u\right)\right)}^2\left(D_{1,2}\right)\left(f\right)\left(t,u\right)D_1\left(f\right)\left(t,u\right)}{3600}+\frac{\left(-15{\left(D_1\left(f\right)\left(t,u\right)\right)}^2\left(D_{2,2}\right)\left(f\right)\left(t,u\right)-30\left(D_{1,1,2}\right)\left(f\right)\left(t,u\right)D_1\left(f\right)\left(t,u\right)-20\left(D_{1,1}\right)\left(f\right)\left(t,u\right)\left(D_{1,2}\right)\left(f\right)\left(t,u\right)-5\left(D_{1,1,1,1}\right)\left(f\right)\left(t,u\right)\right)D_2\left(f\right)\left(t,u\right)}{648000}+\frac{\left(D_{1,2,2}\right)\left(f\right)\left(t,u\right){\left(D_1\left(f\right)\left(t,u\right)\right)}^2}{43200}+\frac{\left(10\left(D_{1,1}\right)\left(f\right)\left(t,u\right)\left(D_{2,2}\right)\left(f\right)\left(t,u\right)+10\left(D_{1,1,1,2}\right)\left(f\right)\left(t,u\right)\right)D_1\left(f\right)\left(t,u\right)}{648000}+\frac{\left(D_{1,1,2}\right)\left(f\right)\left(t,u\right)\left(D_{1,1}\right)\left(f\right)\left(t,u\right)}{64800}+\frac{\left(D_{1,1,1,1,1}\right)\left(f\right)\left(t,u\right)}{648000}]$

结束语

以上。如果有感兴趣朋友，可以进QQ群836204107获取更多Maple相关资源。Maple技术交互研讨群欢迎每一位志同道合朋友！