常微分方程的数值求解及其应用

数值分析设计报告

常微分方程的数值求解及其应用

一、设计目的 .......................................... 2 二、任务简介 .......................................... 2 三、理论基础 .......................................... 3

3.1二阶龙格-库塔方法 ............................. 3 3.2三阶龙格-库塔方法 ............................. 4 3.3四阶龙格-库塔方法 ............................. 5 3.4二阶亚当姆斯格式 .............................. 6 3.5三阶亚当姆斯格式 .............................. 7 3.6四阶亚当姆斯格式 .............................. 8 3.7亚当姆斯预报-校正系统 ......................... 9

四、案例的运算结果 ................................... 10

4.1案例一运算结果 ............................... 10 4.2案例二运算结果 ............................... 11

五、数值分析设计的GUI界面 ........................... 12 六、结果分析 ......................................... 13 七、设计心得 ......................................... 13 八、附录 ............................................. 13

1

一、设计目的

在matlab环境下熟悉的运用计算机编程语言并结合龙格－库塔法、亚当姆斯方法的理论基础对常微分方程组问题进行求解，在运行完程序后以及对运行结果做出各方面的分析和比较。并利用界面将所得结果表示出来。

二、任务简介

用熟悉的计算机语言编程上机完成用二阶龙格－库塔法、三阶龙格－库塔法、四阶龙格－库塔法、二阶亚当姆斯方法、三阶亚当姆斯方法、四阶亚当姆斯方法、亚当姆斯预报--校正系统求解常微分方程组。 设计统一界面的相关算法。

建立模型用数值和解析两种方法解决应用型问题。 需解决的案例如下：

z'1z2案例一：求解微分方程组z'2z3在区间H[0.1,60]上满足条

1233z'3xz33xz22xz19xsinx件：x0.1时，z1z2z31的特解，画出数值解的图像，进行比较。 案例二：放射性废物的处理

美国原子能委员会以往处理浓缩的放射性废料的方法，一直是把它们装入密封的圆桶里，然后扔到水深为90多米的海底。生态学家和科学家们表示担心，怕圆桶下沉到海底时与海底碰撞而发生破裂，从而造成核污染。原子能委员会分辨说这是不可能的。为此工程师们进行了碰撞实验，发现当圆桶下沉速度超过12.2m/s与海底相撞时，圆桶就可能发生碰裂。这样为避免圆桶碰裂，需要计算一下圆桶沉到海底时速度是多少？这时已知圆桶重量为239.46kg，体积为0.2058m3， 海水密度为1035.71 kg /m3。如果圆桶速度小于12.2m/s，就说明这种方法是安全可靠的，否则就要禁止用这种方法来处理放射性废料。假设水的阻力与速度大小成正比例，其正比例常数k=0.6。

（1） 建立解决上述问题的微分方程数学模型。

（2） 用数值和解析两种方法求解微分方程，并回答是否要禁止用这种方法来处理放射性废料。

2

三、理论基础

3.1二阶龙格-库塔方法

yn1ynh[(1)k1k2] k1f(xn,yn)kf(x,yphk)npn12只要成立p1/2，满足这一条件的格式统称为二阶龙格-库塔格式。在本题中

1/2,p1

当1/2,p1时，该式为改进的欧拉格式。

当1,p1/2时，该式为变形的欧拉格式。

二阶龙格库塔-算法流程图：

图3.1 二阶经典的龙格-库塔方法

3

开始 输入x0，y0，h，N n1 x0hx1,f(x0,y0)k1,f(x1,y0hk1)k2y0h/2(k1k2)y1 输出x1,y1 n1nnN? x1x0,y1y0 = = 结束

3.2三阶龙格-库塔方法

类似前面改进的欧拉方法公式的推导，将yn1在(xn,yn)处作Taylor展开，然后再将y(xn1)4xxy(x)yo(h)。于是得nn1n1在处作Taylor展开，只要将两个展开式前四项相同便有

到三阶龙格-库塔公式为：

yn1ynh/6(k14k2k3)kf(x,y)1nn kf(x,yh/2k)n1/2n12k3f(xn1,ynh(k12k2))三阶龙格库塔-算法流程图：

图3.2 三阶经典的龙格-库塔方法

4

开始 输入x0，y0，h，N n1 x0hx1,f(x0,y0)k1,f(x0h/2,y0h/2k1)k2y0h/6(k14k2k3)y1输出x1,y1 n1nf(x0h,y0h(k12*k2))k3 nN? = 结束x1x0,y1y0

3.3四阶龙格-库塔方法

类似前面三阶龙格-库塔的推导方法，如果每步计算四次函数f（x，y）的值，完全类似

5o(h)的四阶龙格-库塔公式，其公式为： 的，可以导出局部截断误差为

yn1ynh/6(k12k22k3k4)kf(x,y)1nn k2f(xn1/2,ynh/2k1)kf(xn1/2,ynh/2k2)3k4f(xn1,ynhk3)四阶龙格库塔-算法流程图：

图3.3四阶经典的龙格-库塔方法

5

开始 输入x0，y0，h，N n1 x0hx1,f(x0,y0)k1,f(x0h/2,y0h/2k1)k2f(x0h/2,y0h/2k2)k3f(x1,y0hk3)k4y0h/6(k12k22k3k4)y1输出x1,y1 nN? n1nx1x0,y1y0 = 结束

3.4二阶亚当姆斯格式

为使计算格式具有二阶精度，须取1/2，这样可推出二阶亚当慕斯格式为：

''yn1ynh/2(3ynyn1)' 类似可推出三阶亚当姆斯格式和四阶亚当姆斯格式。 ynf(xn,yn)'yn1f(xn1,yn1)二阶亚当姆斯方法流程图：

图3.4 二阶亚当姆斯预方法流程图

6

开始 输入x0,y0,h,N 用经典格式（16）提供 y1,y2,y3x0khxk,f(xk,yk)yk'k0,1,2,3 n=2 x1hx2'y1h/2(3y1'y0)y2 f(x1,y1)y1''f(x0,y0)y0输出x2,y2 n1n,x2x1nN? = 结束y2y1,y'k1y,k0,1,2,3'k 

3.5三阶亚当姆斯格式

'''yn1ynh/12(23yn16yn15yn2)' ynf(xn,yn)'yn1f(xn1,yn1)三阶亚当姆斯方法流程图：

图3.5 三阶亚当姆斯方法流程图

7

'''yn1ynh/12(23yn16yn15yn2)'ynf(xn,yn)'yn1f(xn1,yn1)开始 输入x0,y0,h,N 用经典格式（16）提供 y1,y2,y3x0khxk,f(xk,yk)yk'k0,1,2,3 n=3 输出x3,y3  nN? = 结束 n1n,x3x2y3y2,y'k1y,k0,1,2,3'k

3.6四阶亚当姆斯格式

''''yn1ynh/24(55yn59yn137yn29yn3)' ynf(xn,yn)'yn1f(xn1,yn1)四阶亚当姆斯方法流程图：

''''yn1ynh/24(55yn59yn137yn29yn3)'ynf(xn,yn)'yn1f(xn1,yn1)开始 输入x0,y0,h,N 用经典格式（16）提供 y1,y2,y3x0khxk,f(xk,yk)yk'k0,1,2,3 n=4 输出x4,y4 nN? = 结束  n1n,x4x3y4y3,y'k1y,k0,1,2,3'k 图3.6 四阶亚当姆斯方法程图

8

3.7亚当姆斯预报-校正系统 预报：

''''yn1ynh/24(55yn59yn137yn29yn3) 'yn1f(xn1,yn1)校正：

''''yn1ynh/24(9yn119yn5yn1yn2) 'yn1f(xn1,yn1)亚当姆斯预报-校正系统流程图：

图3.7 亚当姆斯预报校正系统流程图

9

开始 输入x0,y0,h,N 用经典格式（16）提供 y1,y2,y3x0khxk,f(xk,yk)yk'k0,1,2,3n=4 x3hx4'''y3h/24(55y359y237y1'9y0)ypf(x4,yp)y'p''y3h/24(9y'p19y35y2y1')y4'f(x4,y4)y4 输出x4,y4 n1n,x4x3nN? = 结束  y4y3,y'k1y,k0,1,2,3'k 四、案例的运算结果

4.1案例一运算结果

二、三、四阶龙格-库塔方法求解方程组所得图像：

2x 1061.510.50-0.5-1-1.5-2-2.50102030405060 图4.1二阶龙格-库塔方法所得图像 2x 1061.510.50-0.5-1-1.5-20102030405060图4.3四阶龙格-库塔方法所得图像 亚当姆斯法求解方程组所得图像：

2.5x 10621.510.50-0.5-1-1.5-20102030405060

图4.4二阶亚当姆斯解得图像：

62x 101.510.50-0.5-1-1.5-20102030405060图4.2三阶龙格-库塔方法所得图像

2x 1061.510.50-0.5-1-1.5-20102030405060

图4.5三阶亚当姆斯解得图像：

10

21.510.50-0.5-1-1.5-2x 10621.510.50-0.5-1-1.5-2x 10601020304050600102030405060

图4.6四阶亚当姆斯解得图像： 图4.7亚当姆斯预报-校正解得图像： 4.2案例二运算结果 4.2.1建立模型

圆桶的位移和速度分别满足下面的微分方程：

d2sds m2mgpgVk （1）

dtdt mdvmgpgVkv （2） dt考虑受到的阻力，这时圆桶的速度应满足如下的微分方程： m4.2.2模型解答 解析法：

根据方程（1）和（2），加上初始条件v|t=0=0，s|t=0=0，以及题设的初始数据。通过Matlab就可以求出圆桶的位移和速度的方程，具体步骤如下：

m=239.46；w=0.2058；g=9.8；p=1035.71；k=0.6；

s=dsolve('m*D2s=m*g-p*g*w-k*Ds','s(0)=0','t') v=dsolve('m*Dv=m*g-p*g*w-k*v','v(0)=0','t')

dvmgpgVkv2 （3） dt 得出位移的方程为：

s（t）= - 171511+429.744t+171511e-0.002505t （4） 速度的方程为：

v（t）= 429.744 - 429.744e-0.002505t （5） 通过方程（4）及s（t）= 90m，利用下面Matlab程序：

solve('90=-171511+429.744*t+171511/exp(0.002505*t)','t')

求出时间t = 12.7868s，再把它代入方程（5），利用Matlab程序：

11

solve('v=429.744-429.744/exp(0.002505*13.002)','v')

求出圆桶的速度为v=13.5504m/s。

显然，此时圆桶的速度已超过12.2m/s，可以得出这种处理废料的方法不合理，美国原子能委员会已经禁止用这种方法来处理放射性废料。

根据方程（3），加上初始条件v| t=0 = 0，可利用下面的Matlab程序：

v=dsolve('m*Dv=m*g-p*g*w-k*v^2','v(0)=0','t') 求出圆桶的速度：v（t）=

FFkTanht。 km 若把题设中F，k（仍设为0.6的话）和m的值代入上式，可得圆桶的速度为： v（t）=20.7303Tanh（0.0519426t）

这时若该速度要小于12.2m/s，那么利用Matlab可得圆桶的运动时间就不能超过11.4940s，位移不能超过70.44m。对应的Matlab源程序如下：

solve('12.2=429.744-429.744/exp(0.002505*t)','t') （求时间）

solve('s=-171511+429.744*11.4940+171511/exp(0.002505*11.4940)','s')

（求位移）

从这个模型，我们也可以得出原来处理核废料的方法是不合理的。 数值法：

求解出函数s ( t ) 后，求出圆桶下落到水深为9 0 米海底的时间，即求满足s( t ) = 9 0 的时间t， 然后求出在时间t 的速度即可以得到问题的解答。这里的速度关系由直接求函数s ( t ) 对t 的导数得到。

画出方程s( t ) = 9 0 的函数图形找根t, 可以发现根t在0,20内，继续画出s(t)在0,20内的局部图形，不断代入不同的t的值，可以发现根t在14附近，再将t14带入v(t)的表达式中，得出v(t)。结果表明当下降深度为90米的时候，速度接近于13.5米每秒，大于12.2米每秒。因此，计算结果说明美国原子能委员会以往处理浓缩的放射性废料方法是不安全的。

五、数值分析设计的GUI界面

图5.1从GUI界面所得结果图像

12

六、结果分析

由以上的各种方法与精确解的比较图可以看出在计算精度上，四阶龙格-库塔方法和亚当姆斯预报-校正系统的误差最小，其次是三阶龙格-库塔方法和四阶亚当姆斯方法，再次之是三阶亚当姆斯方法，二阶龙格-库塔方法和二阶亚当姆斯方法误差则比较大，所以四阶龙格-库塔方法和亚当姆斯预报-校正系统得到最佳的精度。而在计算量上面，相应地，很明显的四阶龙格-库塔方法和亚当姆斯预报-校正系统也是最大，二阶龙格-库塔方法和二阶亚当姆斯方法计算量最小。这样的结果，说明了运用以上三种方法时，其计算量的多少与精度的大小成正比。我们在实际运用与操作中，可以根据实际情况，选择这上述方法中的其中一种最适合的，追求精度的话，可以使用四阶龙格-库塔方法和亚当姆斯预报-校正系统；而二阶龙格-库塔方法和二阶亚当姆斯方法，在精度上和计算量上都不错，能够满足一般情况；

七、设计心得

这次实验花了较多的时间，先是选择方程方面遇到了困难：有些方程在运行过程中会不断出错，错误的原因有很多，有的是在输入格式上有错，有的是在步长和初值的取法上出错。在试了十多个方程之后终于找到一个自己满意的。接下来的运行中也出现了一些错误。有的方程即使求的出来结果与精确解的误差是很大的，显然是错误的。还有就是刚开始犯了一个低级错误：把原函数直接代入各种格式求解这样当然是不正确的。取初值对函数图像的影响很大，如果初值取的不合适误差也会非常大，所以要合理的选择初值，在几种求解方法中只有龙格—库塔的求解与精确值最为接近几乎是重合的。通过的后面的不同步长求解的比较发现步长越小与精确值越接近。总之在做这次实验中对学习数学有了更多的体会，那就是一定要认真，不能丝毫马虎，一个小小的错误就会导致所有的数据不能运行。其做的主要是亚当姆斯算法，做的主要是龙格-库塔算法，至于界面管理和案例二我们俩个都有参与，共同完成。这也更好的锻炼了我们的团队合作的能力。

八、附录

二阶龙格-库塔方法：

function [x,y]=marunge2s(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2);

y=zeros(length(y0),length(x)); y(:,1)=y0(:); for n=1:(length(x)-1)

13

k1=feval(dyfun,x(n),y(:,n)); k2=feval(dyfun,x(n+1),y(:,n)+h*k1); y(:,n+1)=y(:,n)+h*(1/2*k1+1/2*k2); end

function ff=mafun(x,y) ff(1)=y(2); ff(2)=y(3);

ff(3)=y(3)/x-3*y(2)/(x^2)+2*y(1)/(x^3)+9*x^3*sin(x); ff=ff(:);

三阶龙格-库塔方法：

function [x,y]=marunge3s(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2);

y=zeros(length(y0),length(x)); y(:,1)=y0(:); for n=1:(length(x)-1) k1=feval(dyfun,x(n),y(:,n));

k2=feval(dyfun,x(n)+h/2,y(:,n)+h/2*k1); k3=feval(dyfun,x(n+1),y(:,n)+h*(-k1+2*k2)); y(:,n+1)=y(:,n)+(h/6).*(k1+4*k2+k3); end

function ff=mafun(x,y) ff(1)=y(2); ff(2)=y(3);

ff(3)=y(3)/x-3*y(2)/(x^2)+2*y(1)/(x^3)+9*x^3*sin(x); ff=ff(:);

四阶龙格-库塔方法：

function [x,y]=marunge4s(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2);

y=zeros(length(y0),length(x)); y(:,1)=y0(:); for n=1:(length(x)-1) k1=feval(dyfun,x(n),y(:,n));

k2=feval(dyfun,x(n)+h/2,y(:,n)+h/2*k1); k3=feval(dyfun,x(n)+h/2,y(:,n)+h/2*k2); k4=feval(dyfun,x(n+1),y(:,n)+h*k3); y(:,n+1)=y(:,n)+(h/6)*(k1+2*k2+2*k3+k4);

14

end

function ff=mafun(x,y) ff(1)=y(2); ff(2)=y(3);

ff(3)=y(3)/x-3*y(2)/(x^2)+2*y(1)/(x^3)+9*x^3*sin(x); ff=ff(:);

二阶亚当姆斯方法：

function [x,y]=yadangmusi2s(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2);

y=zeros(length(y0),length(x)); y(:,1)=y0(:); for n=2:(length(x)-1) k1=feval(dyfun,x(n),y(:,n)); k2=feval(dyfun,x(n-1),y(:,n-1)); y(:,n+1)=y(:,n)+(h/2)*(3*k1-k2); end

function ff=mafun(x,y) ff(1)=y(2); ff(2)=y(3);

ff(3)=y(3)/x-3*y(2)/(x^2)+2*y(1)/(x^3)+9*x^3*sin(x); ff=ff(:);

三阶亚当姆斯方法：

function [x,y]=yadangmusi3s(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2);

y=zeros(length(y0),length(x)); y(:,1)=y0(:); for n=3:(length(x)-1) k1=feval(dyfun,x(n),y(:,n)); k2=feval(dyfun,x(n-1),y(:,n-1)); k3=feval(dyfun,x(n-2),y(:,n-2));

y(:,n+1)=y(:,n)+(h/12)*(23*k1-16*k2+5*k3); end

function ff=mafun(x,y) ff(1)=y(2); ff(2)=y(3);

ff(3)=y(3)/x-3*y(2)/(x^2)+2*y(1)/(x^3)+9*x^3*sin(x);

15

ff=ff(:);

四阶亚当姆斯方法：

function [x,y]=yadangmusi4s(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2);

y=zeros(length(y0),length(x)); y(:,1)=y0(:); for n=4:(length(x)-1) k1=feval(dyfun,x(n),y(:,n)); k2=feval(dyfun,x(n-1),y(:,n-1)); k3=feval(dyfun,x(n-2),y(:,n-2)); k4=feval(dyfun,x(n-3),y(:,n-3));

y(:,n+1)=y(:,n)+(h/24)*(55*k1-59*k2+37*k3-9*k4); end

function ff=mafun(x,y) ff(1)=y(2); ff(2)=y(3);

ff(3)=y(3)/x-3*y(2)/(x^2)+2*y(1)/(x^3)+9*x^3*sin(x); ff=ff(:);

亚当姆斯预报-校正系统：

function [x,y]=yadangmusi5s(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2);

y=zeros(length(y0),length(x)); y(:,1)=y0(:); for n=4:(length(x)-1) k1=feval(dyfun,x(n),y(:,n)); k2=feval(dyfun,x(n-1),y(:,n-1)); k3=feval(dyfun,x(n-2),y(:,n-2)); k4=feval(dyfun,x(n-3),y(:,n-3));

k0=y(:,n)+(h/24)*(55*k1-59*k2+37*k3-9*k4); k=feval(dyfun,x(n+1),k0);

y(:,n+1)=y(:,n)+(h/24)*(9*k+19*k1-5*k2+k3); end

function ff=mafun(x,y) ff(1)=y(2); ff(2)=y(3);

ff(3)=y(3)/x-3*y(2)/(x^2)+2*y(1)/(x^3)+9*x^3*sin(x);

16

ff=ff(:);

GUI界面aniu:

function pushbutton4_Callback(hObject, eventdata, handles) [x,y]=marunge2s(@mafun,[0.1 60],[1 1 1],0.4); plot(x,y(1,:),x,y(2,:),x,y(3,:),'r');

function pushbutton5_Callback(hObject, eventdata, handles) [x,y]=marunge3s(@mafun,[0.1 60],[1 1 1],0.4); plot(x,y(1,:),x,y(2,:),x,y(3,:),'r');

function pushbutton6_Callback(hObject, eventdata, handles) [x,y]=marunge4s(@mafun,[0.1 60],[1 1 1],0.4); plot(x,y(1,:),x,y(2,:),x,y(3,:),'r');

function pushbutton7_Callback(hObject, eventdata, handles) [x,y]=yadangmusi2s(@mafun,[0.1 60],[1 1 1],0.4); plot(x,y(1,:),x,y(2,:),x,y(3,:),'r');

function pushbutton8_Callback(hObject, eventdata, handles) [x,y]=yadangmusi3s(@mafun,[0.1 60],[1 1 1],0.4); plot(x,y(1,:),x,y(2,:),x,y(3,:),'r');

function pushbutton9_Callback(hObject, eventdata, handles) [x,y]=yadangmusi4s(@mafun,[0.1 60],[1 1 1],0.4); plot(x,y(1,:),x,y(2,:),x,y(3,:),'r');

function pushbutton10_Callback(hObject, eventdata, handles) [x,y]=yadangmusi5s(@mafun,[0.1 60],[1 1 1],0.4); plot(x,y(1,:),x,y(2,:),x,y(3,:),'r');

17