数值计算方法后退欧拉法,龙格库塔法(3阶,4阶方法).docVIP

数值计算方法实验课测验作业——微分方程求解的后退欧拉法,龙格库塔法（三阶，四阶方法） 日期： 2011-06-17 一、 实验目的 1. 学习matlab的使用方法。 2. 掌握常微分方程的几种数值解法：后退欧拉法，龙格库塔法三阶方法，龙格库塔法四阶方法。 3. 比较各方法的数值解及误差，了解各方法的优缺点。 二、 实验题目 给定的初值问题 y′=-y+x+2,0≤x≤1 y(0)=-1, 取精确解 y(x)=exp(-x)+x 按 (1) 后退欧拉法，步长h=0.003, h=0.1; (2) 龙格库塔法三阶方法，步长h=0.1; (3) 龙格库塔法四阶方法，步长h=0.1; 求在节点xk=1+0.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+2 end ②精确解计算(ouleij_q) function z=ouleij_q(x) z=exp(-x)+x end 2.对于龙哥库塔三阶： 利用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=-1; x=1; for i=1:20 k1=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.2084 z =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.5152 对于龙哥库塔三阶： ①若 clear all h=0.1; y=-1; x=1; for i=1:25; k1=oulei_wf(x,y);