数值计算方法后退欧拉法龙格库塔法三阶四阶方法

数值计算方法实验课测验作业——微分方程求解的后退欧拉法,龙格库塔法(三阶，四阶方法)日期：2011-06-17、实验目的1.2.学习matlab的使用方法掌握常微分方程的几种数值解法：后退欧拉法，龙格库塔法三阶方法，龙格库 塔法四阶方法3.比较各方法的数值解及误差，了解各方法的优缺点 、实验题目给定的初值问题y' =-y+x+2 , Ow x < 1 y(o)=-1,取精确解 y(x)=exp(-x)+ x按 ⑴后退欧拉法，步长 h=0.003, h=0.1;(2) 龙格库塔法三阶方法，步长 h=0.1;(3) 龙格库塔法四阶方法，步长 h=0.1;求在节点Xk=1+O.1k (k=1,2,3 10)处的数值解及误差比较各方法的优缺点三、实验原理1. 对于后退欧拉法：利用 Yn+1 = Yn + 1/2*K1 + 1/2*K2 ① n = 1, 2, 3 ……K1 = hf (Xn, Yn ) ②K2 = hf (Xn + h, Yn + K1 ) ③ 三式可以完成计算需要将微分方程表达式和精度计算表达式作为两个函数保存在 m文件里并在程序中调用：① 微分方程(oulei_ wf)function z=oulei_wf(x,y)z=-y+x+2end② 精确解计算(ouleij_q)function z=ouleij_q(x)z=exp(-x)+xend2. 对于龙哥库塔三阶：禾U用 Yn+1 = Yn + h/6(K1 + 4K2 + K3 ) ①K1 :=f( Xn, Yn)②K2 ==f ( Xn + 1/2*h, Yn +( h/2)*K1)③K3= f( Xn + h, Yn - K1+2*h*k2 )④四式可以完成计算3.对于龙哥库塔四阶：利用 Yn+1 = Yn + 1/6(K1 + 2K2 + 2K3 +K4)①K1 ==hf ( Xn, Yn)②K2 ==hf ( Xn + 1/2*h, Yn + 1/2*K1)③K3 ==hf ( Xn + 1/2*h, Yn + 1/2*K2)④K4 ==hf ( Xn + h, Yn + K3)⑤四式可以完成计算四、 实验内容由上述实验原理叙述的后退欧拉法， 龙哥库塔三阶，龙哥库塔 四阶几种常微分方程数值解法分别对已给定的初值问题进行 求解，比较各方法的数值解及误差，了解各方法的优缺点。

五、 实验结果1. 对于后退欧拉法：①若>>h=0・1;y=-i;x=1;for i=1:20k1=h*oulei_wf(x,y); k2=h*oulei_wf(x+h,y+k1); y=y+0 ・5*k1+0 ・5*k2 x=x+h;z=ouleij_q(x)t=y-z end z = 4,z = 3.7000,y = -0.6150 z = 1.4329, z =1.4329, t = -2.0479 z =3.7150, z =3.4435, y = -0.2571 z =1.5012, z =1.5012, t =-1.7583 z =3.4571, z = 3.2114, y =0.0763 z =1.5725, z =1.5725, t =-1.4962 z =3.2237, z = 3.0013, y =0.3876 z =1.6466, z =1.6466, t =-1.2590 z =3.0124,z =2.8112,y =0.6788 z =1.7231,z =1.7231,t =-1.0444 z =2.8212,z = 2.6391,y =0.9518 z =1.8019,z =1.8019,t =-0.8501 z = 2.6482,z =2.4834,y =1.2084z =1 ・8827,z =1.8827, t =-0.6743 z =2.4916,z =2.3425, y =1.4501 z =1.9653,z =1.9653, t =-0.51522. 对于龙哥库塔三阶：①若 >> clear allh=0.1;y=-1;x=1;for i=1:25;k1=oulei_wf(x,y); k2=oulei_wf(x+h/2,y+k1*(h/2)); k3=oulei_wf(x+h,y-k1*h+k2*(2*h)); y=y+(k1+4*k2+k3)*(h/6);x=x+h;z=ouleij_q(x);t=abs(y-z);A=[x y z t]endz = 4,z = 3.8500,z = 3.7300,z =1.4329A = 1.1000 -0.6145 1.4329 2.0474z =3.7145z=3.5788z = 3.4702z =1.5012 A = 1.2000 -0.2562 1.5012 1.7574z =3.4562z =3.3334z =3.2351z =1.5725 A = 1.3000 0.0776 1.5725 1.4950z = 3.2224z =3.1113z =3.0224z =1.6466 A = 1.4000 0.3891 1.6466 1.2575z = 3.0109z = 2.9104z =2.8299z =1.7231 A =1.5000 0.6804 1.7231 1.0427z =2.8196z =2.7286z =2.6558z =1.8019 A = 1.6000 0.9536 1.8019 0.8483z =2.6464z =2.5641z =2.4982z =1.8827 A = 1.7000 1.2103 1.8827 0.6724z =2.4897z =2.4152z =2.3556z =1.9653 A =1.8000 1.4521 1.9653 0.5132z =2.3479z =.2805z =2.2266z =2.0496 A =1.9000 1.6803 2.0496 0.36923. 龙哥库塔四阶： >> clear all h=0.1;y=-1;x=1;for i=1:25;k1=oulei_wf(x,y);k2=oulei_wf (x+h/2,y+k1*(h/2));k3=oulei_wf (x+h/2,y+k2*(h/2));k4=oulei_wf (x+h,y+h*k3);y=y+(k1+2*k2+2*k3+k4)*(h/6);x=x+h;z=ouleij_q(x);t=abs(y-z);A=[x y z t]endz =4z =3.8500z =3.8575z =3.7142z =1.4329A =1.1000 -0.6145 1.4329 2.0474z =3.7145z =3.5788z =3.5856z =3.4560z =1.5012A =1.2000 -0.2562 1.5012 1.7574z =3.4562z =3.3334z =3.3395z =3.2222z =1.5725A =1.3000 0.0775 1.5725 1.4950z =3.2225z =3.1113z =3.1169z =3.0108z =1.6466A =1.4000 0.3890 1.6466 1.2576z =3.0110z =2.9104z =2.9154z =2.8194z =1.7231A =1.5000 0.6804 1.7231 1.0427z =2.8196z =2.7286z =2.7332z =2.6463z =1.8019A =1.6000 0.9536 1.8019 0.8483z =2.6464z =2.5641z =2.5682z =2.4896z =1.8827A =1.7000 1.2102 1.8827 0.6724z =2.4898z =2.4153z =2.4190z =2.3479z =1.9653A =1.8000 1.4520 1.9653 0.5133z =2.3480z =2.2806z =2.2840z =2.2196z =2.0496A =1.9000 1.6803 2.0496 0.3693z =2.2197z =2.1587z =2.1618z =2.1035z =2.1353A =2.0000 1.8964 2.1353 0.2390z =2.1036z =2.0485z =2.0 512z =1.9985z =2.2225A =2.1000 2.1014 2.2225 0.1211z =1.9986z =1.9487z =1.9512z =1.9035z =2.3108A =2.2000 2.2964 2.3108 0.0144六、实验结果分析1. 对于欧拉法，步长越小，精度越高，而产生的误差越小。

总 体来说，欧拉法的优点是形式简单，计算方便，缺点是总的运 算精度比较低而且随着x的增大，误差值也越来越大根据 欧拉公式的截断误差计算，欧拉法是一阶方法2. 对于改进欧拉法，其基本特征与欧拉法相似，也是步长越小，精度越咼，误差越小优点是精度相对欧拉法来说较咼， 截断误差为0(h^3)，缺点是比欧拉法计算量大根据欧拉改进法 公式的截断误差计算，欧拉改进法是二阶方法3. 对于龙格一库塔法，优点是精度更咼，同样的步长下精度比 欧拉法高的多，误差更小，截断误差为0(h^5)，故龙格一库 塔法是四阶方法缺点是每步都要计算四次微分值4. 综上，选取方法时，可综合考虑精度要求和复杂度控制要求 等实际需要，从而选择适当的方法求解。