龙格库塔法解微分方程Maple

这是第一次写文章，有点胆怯，发一篇两年前写的关于如何用Maple基于龙格库塔法求解一类简单常微分方程。由于Maple内置dsolve命令已经十分强势，这里更多是体验方法本身。理论部分随手找一本ode书籍就OK，现在我们直接上代码：

#Janeasefor do it in 2017/12/29 with Maple2017~
restart:
with(plots):

#function content
four_stage_Runge_Kutta_dsolve:=proc(func::procedure,dt,n::integer,u0)
description"dsolve ordinary differential equation like diff(u(t),t)=func(t,u)";
local a,b,c,k,K,K_expr,tl,ul,i;
a:=[0,1/3,2/3,1];b:=[0,[1/3],[-1/3,1],[1,-1,1]];c:=[1/8,3/8,3/8,1/8];#Kutta methed
#a:=[0,1/2,1/2,1];b:=[0,[1/2],[0,1/2],[0,0,1]];c:=[1/6,2/6,2/6,1/6];#classic Runge_Kutta methed
tl:=[seq(ha*dt,ha=1..n)];ul:=tl;
k[1]:=func(t,u);
for i in [2,3,4] do
    k[i]:=func(t+a[i]*dt,u+dt*sum(b[i][ha]*k[ha],ha=1..i-1));
end do;
K_expr:=simplify(sum(c[ha]*k[ha],ha=1..4));
K:=(kt,ku)->eval(K_expr,[t=kt,u=ku]);
ul[1]:=evalf(u0+dt*K(0,u0),15);
for i in [$2..n] do
    ul[i]:=evalf(ul[i-1]+dt*K(tl[i-1],ul[i-1]),15);
end do;
tl,ul;
end proc:

#example and pictures by contrast
func:=(t,u)->-2*u;dt:=0.1;n:=10;u0:=1;
ans:=four_stage_Runge_Kutta_dsolve(func,dt,n,u0):
display(plot(ans),plot(u0*exp(-2*x),x=0..1,color=blue),legend=["Runge_kutta","real"],title=typeset("Runge_Kutta_method's test on",diff(u(t),t)=func(t,u)));
errorvalue:=map(x->u0*exp(-2*x),ans[1])-ans[2];

我这本参考书系数提到经典龙格库塔法与龙格库塔法两种对应系数，这里我都写进程序里面，无非一个用备注#注释掉，必要时注释掉就可以。还有在ode输入方面我采用函数格式而非方程格式（最好是方程格式，这里由于更多体验方法所以侧重个人程序书写便利而采用函数格式）。我们以下面方程为例：

结果对比图如下：

可以发现，二者重合得很好（当然除了理论优越性外我的测试函数找的好，如果本来就是发散函数注定是不会有结果的~） 

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