微分方程求解数值方法.doc

重庆三峡学院数值计算方法实验课测验微分方程求解的后退欧拉法（后退梯形法）、龙格库塔法（三阶、四阶公式） 院 系 数学与统计学院 专 业 信息与计算科学 姓 名 庞欢 年 级 2009级 学 号 200906014132 指导教师 刘学飞 2011年6月17日数值计算方法实验测试作业 微分方程的后退欧拉法、龙格库塔法（三阶、四阶方法） 一、实验目的1. 学习matlab的使用方法。2.掌握常微分方程的几种数值解法：后退欧拉法，龙格库塔法三阶方法，龙格库塔法四阶方法。3. 比较各方法的数值解及误差，了解各方法的优缺点。二、实验题目给定的初值问题： 精确解计算公式： 按 (1)后退欧拉法，步长h=0.025, h=0.1; (2)三阶标准龙格库塔法，步长h=0.02; (3)四阶标准龙格库塔法，步长h=0.02;求在各节点处的数值解及误差，并比较各方法的优缺点。三、实验原理1.对于后退欧拉法：利用Yn+1 = Yn + 1/2*K1 + 1/2*K2 n = 1, 2, 3K1 = hf（Xn, Yn） K2 = hf（Xn + h, Yn + K1） 三式可以完成计算需要将微分方程表达式和精度计算表达式作为两个函数保存在m文件里并在程序中调用： 微分方程（pang_ huan）function z=pang_huan(x,y)z=x*x+x-y;end 精确解计算(panghuanj_q)function z=panghuanj_q(x)z=-exp(-x)+x*x-x+1end 2、对于龙格库塔三阶解法：利用Yn+1 = Yn + 1/2*K1 + 1/2*K2 n = 1, 2, 3K1 = hf（Xn, Yn） K2 = hf（Xn + h, Yn + K1） 三式可以完成计算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.025; y=0;x=0;for i=1:40;y=y+h*pang_huan(x,y)x=x+hvpa(y)z=panghuanj_q(x);t=abs(y-z);x0=x;plot(x0,z,*,x0,y)A=x y z tend% x0=1:0.025:1% z0=panghuanj_q(x0)运行结果：A = 0.0250 0 0.0003 0.0003A = 0.0500 0.0006 0.0013 0.0006A = 0.0750 0.0019 0.0029 0.0009A = 0.1000 0.0039 0.0052 0.0013A= 0.1250 0.0066 0.0081 0.0016A = 0.1500 0.0099 0.0118 0.0019A = 0.1750 0.0140 0.0162 0.0022A = 0.2000 0.0188 0.0213 0.0025A = 0.2250 0.0243 0.0271 0.0028A = 0.2500 0.0306 0.0337 0.0031A = 0.2750 0.0376 0.0411 0.0034A = 0.3000 0.0455 0.0492 0.0037A = 0.3250 0.0541 0.0581 0.0040A = 0.3500 0.0635 0.0678 0.0043A = 0.3750 0.0737 0.0783 0.0046A= 0.4000 0.0848 0.0897 0.0049A = 0.4250 0.0966 0.1019 0.0052A = 0.4500 0.1094 0.1149 0.0055A = 0.4750 0.1229 0.1287 0.0058A = 0.5000 0.1374 0.1435 0.0062A = 0.5250 0.1527 0.1591 0.0064A = 0.5500 0.1689 0.1756 0.0067A = 0.5750 0.1860 0.1929 0.0069A = 0.6000 0.2040 0.2112 0.0072A = 0.6250 0.2229 0.2304 0.0075A = 0.6500 0.2427 0.2505 0.0078A = 0.6750 0.2634 0.2715 0.0080A = 0.7000 0.2851 0.2934 0.0083A = 0.7250 0.3077 0.3163 0.0086A = 0.7500 0.3313 0.3401 0.0088A = 0.7750 0.3558 0.3649 0.0091A = 0.8000 0.3813 0.3907 0.0093A = 0.8250 0.4078 0.4174 0.0096A = 0.8500 0.4352 0.4451 0.0098A = 0.8750 0.4637 0.4738 0.0101A = 0.9000 0.4931 0.5034 0.0103A = 0.9250 0.5235 0.5341 0.0106A = 0.9500 0.5550 0.5658 0.0108A = 0.9750 0.5874 0.5984 0.0110A = 1.0000 0.6208 0.6321 0.0113结果分析：由运行可看出，后退欧拉法与最基础的欧拉法是相似的。但缺点是x的增大误差也越大，且运算繁杂。2、对于龙格库塔的三阶解法：h=0.02; y=0; x=0;for i=1:50;k1=h*pang_huan(x,y);k2=h* pang_huan (x+h/2,y+k1/2);k3=h* pang_huan(x+h/2,y+k2/2);k4=h*pang_huan(x+h,y+k3);y=y+(k1+2*k2+2*k3+k4)/6;x=x+h;z=panghuanj_q(x);t=abs(y-z);A=x y z tEnd运行结果：A = 0.0200 0.0002 0.0002 0.0000A = 0.0400 0.0008 0.0008 0.0000A = 0.0600 0.0018 0.0018 0.0000A = 0.0800 0.0033 0.0033 0.0000A = 0.1000 0.0052 0.0052 0.0000A = 0.1200 0.0075 0.0075 0.0000A = 0.1400 0.0102 0.0102 0.0000A = 0.1600 0.0135 0.0135 0.0000A = 0.1800 0.0171 0.0171 0.0000A = 0.2000 0.0213 0.0213 0.0000A = 0.2200 0.0259 0.0259 0.0000A = 0.2400 0.0310 0.0310 0.0000A = 0.2600 0.0365 0.0365 0.0000A = 0.2800 0.0426 0.0426 0.0000A = 0.3000 0.0492 0.0492 0.0000A = 0.3200 0.0563 0.0563 0.0000A = 0.3400 0.0638 0.0638 0.0000A = 0.3600 0.0719 0.0719 0.0000A = 0.3800 0.0805 0.0805 0.0000A = 0.4000 0.0897 0.0897 0.0000A = 0.4200 0.0994 0.0994 0.0000A = 0.4400 0.1096 0.1096 0.0000A = 0.4600 0.1203 0.1203 0.0000A = 0.4800 0.1316 0.1316 0.0000A = 0.5000 0.1435 0.1435 0.0000A = 0.5200 0.1559 0.1559 0.0000A = 0.5400 0.1689 0.1689 0.0000A = 0.5600 0.1824 0.1824 0.0000A = 0.5800 0.1965 0.1965 0.0000A = 0.6000 0.2112 0.2112 0.0000A = 0.6200 0.2265 0.2265 0.0000A = 0.6400 0.2423 0.2423 0.0000A = 0.6600 0.2587 0.2587 0.0000A = 0.6800 0.2758 0.2758 0.0000A = 0.7000 0.2934 0.2934 0.0000A = 0.7200 0.3116 0.3116 0.0000A = 0.7400 0.3305 0.3305 0.0000A = 0.7600 0.3499 0.3499 0.0000A = 0.7800 0.3700 0.3700 0.0000A = 0.8000 0.3907 0.3907 0.0000A = 0.8200 0.4120 0.4120 0.0000A = 0.8400 0.4339 0.4339 0.0000A = 0.8600 0.4564 0.4564 0.0000A = 0.8800 0.4796 0.4796 0.0000A = 0.9000 0.5034 0.5034 0.0000A = 0.9200 0.5279 0.5279 0.0000A = 0.9400 0.5530 0.5530 0.0000A = 0.9600 0.5787 0.5787 0.0000A = 0.9800 0.6051 0.6051 0.0000A = 1.0000 0.6321 0.6321 0.0000结果分析：上面是三阶标准龙格库塔法运行的结果，其精度比后退欧拉法高，误差较小。3、对于龙格库塔的四阶解法：h=0.02; y=0; x=0;for i=1:50;k1=pang_huan(x,y);k2=pang_huan (x+h/2,y+k1*(h/2);k3=pang_huan(x+h/2,y+k2*(h/2);k4=pang_huan(x+h,y+h*k3);y=y+(k1+2*k2+2*k3+k4)*(h/6);x=x+h;z=panghuanj_q(x);t=abs(y-z);A=x y z tEnd运行结果：A = 0.0200 0.0002 0.0002 0.0000A = 0.0400 0.0008 0.0008 0.0000A = 0.0600 0.0018 0.0018 0.0000A = 0.0800 0.0033 0.0033 0.0000A = 0.1000 0.0052 0.0052 0.0000A = 0.1200 0.0075 0.0075 0.0000A = 0.1400 0.0102 0.0102 0.0000A = 0.1600 0.0135 0.0135 0.0000A = 0.1800 0.0171 0.0171 0.0000A = 0.2000 0.0213 0.0213 0.0000A = 0.2200 0.0259 0.0259 0.0000A = 0.2400 0.0310 0.0310 0.0000A = 0.2600 0.0365 0.0365 0.0000A = 0.2800 0.0426 0.0426 0.0000A = 0.3000 0.0492 0.0492 0.0000A = 0.3200 0.0563 0.0563 0.0000A = 0.3400 0.0638 0.0638 0.0000A = 0.3600 0.0719 0.0719 0.0000A = 0.3800 0.0805 0.0805 0.0000A = 0.4000 0.0897 0.0897 0.0000A = 0.4200 0.0994 0.0994 0.0000A = 0.4400 0.1096 0.1096 0.0000A = 0.4600 0.1203 0.1203 0.0000A = 0.4800 0.1316 0.1316 0.0000A = 0.5000 0.1435 0.1435 0.0000A = 0.5200 0.1559 0.1559 0.0000A = 0.5400 0.1689 0.1689 0.0000A = 0.5600 0.1824 0.1824 0.0000A = 0.5800 0.1965 0.1965 0.0000A = 0.6000 0.2112 0.2112 0.0000A = 0.6200 0.2265 0.2265 0.0000A = 0.6400 0.2423 0.2423 0.0000A = 0.6600 0.2587 0.2587 0.0000A = 0.6800 0.2758 0.2758 0.0000A = 0.7000 0.2934 0.2934 0.0000A = 0.7200 0.3116 0.3116 0.0000A = 0.7400 0.3305 0.3305 0.0000A = 0.7600 0.3499 0.3499 0.0000A = 0.7800 0.3700 0.3700 0.0000A = 0.8000 0.3907 0.3907 0.0000A = 0.8200 0.4120 0.4120 0.0000A = 0