1、常微分方程课程实习福建农林大学计算机与信息学院（数学类课程）课程实习报告课程名称：常微分方程课程实习实习题目：常微分方程数值求解问题的实习姓 名：系：应用数学专 业：数学与应用数学年 级：2008学 号：指导教师：职 称：讲师2009 年 12 月 1 日福建农林大学计算机与信息学院数学类课程实习报告结果评定评语：成绩：指导教师签字：评定日期：目 录1. 实习的目的和任务 12. 实习要求 13. 实习地点 14. 主要仪器设备 15. 实习内容 15.1 用不同格式对同一个初值问题的数值求解及其分析.15.1.1求精确解 15.1.2用欧拉法求解 35.1.3用改进欧拉法求解 55.1.4用
2、4级4阶龙格库塔法求解 75.1.5 问题讨论与分析95.2 一个算法不同不长求解同一个初值问题及其分析.136. 结束语 22参考文献 22常微分方程课程实习1. 实习的目的和任务目的：通过课程实习能够应用MATLAB软来计算微分方程（组）的数值解；了解常微分方程数值解。任务：通过具体的问题，利用MATLAB软件来计算问题的结果，分析问题的结论。2. 实习要求能够从案例的自然语言描述中，抽象出其中的数学模型；能够熟练应用所学的数值解计算方法；能够熟练使用MATLAB软件；对常微分方程数值解有所认识，包括对不同算法有所认识和对步长有所认识。3. 实习地点数学实验室4. 主要仪器设备计算机、Mi
3、crosoft Windows 7Matlab R2009a5. 实习内容5.1 用欧拉方法，改进欧拉方法，4阶龙格库塔方法分别求下面微分方程的初值dy/dx=y-x2+1 y(0)=0.5 x0，25.1 .1求精确解首先可以求得其精确解为：y=2*x-exp(x)/2 +x2+15.1.1 程序代码： x=0:0.1:2;y=2*x-exp(x)/2+x.2+1plot(x,y,r*-)Data=x,yy = Columns 1 through 17 0.5000 0.6574 0.8293 1.0151 1.2141 1.4256 1.6489 1.8831 2.1272 2.3802
4、2.6409 2.9079 3.1799 3.4554 3.7324 4.0092 4.2835 Columns 18 through 21 4.5530 4.8152 5.0671 5.3055Data = 0 0.5000 0.1000 0.6574 0.2000 0.8293 0.3000 1.0151 0.4000 1.2141 0.5000 1.4256 0.6000 1.6489 0.7000 1.8831 0.8000 2.1272 0.9000 2.3802 1.0000 2.6409 1.1000 2.9079 1.2000 3.1799 1.3000 3.4554 1.40
5、00 3.7324 1.5000 4.0092 1.6000 4.2835 1.7000 4.5530 1.8000 4.8152 1.9000 5.06712.0000 5.30555.1.2 用欧拉法求解程序如下：建立函数文件Euler.mfunction x,y=Euler(fun,x_span,y0,h)x=x_span(1):h:x_span(2);y(1)=y0;for n=1:length(x)-1 y(n+1)=y(n)+h*feval(fun,x(n),y(n);endx=x;y=y;在MATLAB输入以下程序： clear all fun=inline(y-x2+1); x
6、,y=Euler(fun,0,2,0.5,0.1); x,y plot(x,y,b*-)结果及其图象：ans = 0 0.5000 0.1000 0.6500 0.2000 0.8140 0.3000 0.9914 0.4000 1.1815 0.5000 1.3837 0.6000 1.5971 0.7000 1.8208 0.8000 2.0538 0.9000 2.2952 1.0000 2.5438 1.1000 2.7981 1.2000 3.0569 1.3000 3.3186 1.4000 3.5815 1.5000 3.8437 1.6000 4.1030 1.7000 4.3
7、573 1.8000 4.6040 1.9000 4.84052.0000 5.06355.1.3用改进欧拉法求解：程序如下：建立函数文件gjEuler.mfunction x,y= gjEuler (fun,x_span,y0,h)x=x_span(1):h:x_span(2);y(1)=y0;for n=1:length(x)-1 k1=feval(fun,x(n),y(n); y(n+1)=y(n)+h*k1; k2=feval(fun,x(n+1),y(n+1); y(n+1)=y(n)+h*(k1+k2)/2;endx=x;y=y;在MATLAB输入以下程序： clear all f
8、un=inline( y-x2+1); x,y=gjEuler(fun,0,2,0.5,0.1); x,y plot(x,y,m+-)结果及其图象： 0 0.5000 0.1000 0.6570 0.2000 0.8284 0.3000 1.0137 0.4000 1.2122 0.5000 1.4232 0.6000 1.6459 0.7000 1.8794 0.8000 2.1228 0.9000 2.3750 1.0000 2.6348 1.1000 2.9010 1.2000 3.1720 1.3000 3.4464 1.4000 3.7223 1.5000 3.9978 1.6000
9、 4.2708 1.7000 4.5390 1.8000 4.7996 1.9000 5.0499 2.0000 5.28665.1.4 用4阶龙格库塔求解程序如下：建立函数文件Runge_Kutta.mfunction x,y= Runge_Kutta (fun,x_span,y0,h)x=x_span(1):h:x_span(2);y(1)=y0;for n=1:length(x)-1 k1=feval(fun,x(n),y(n); k2=feval(fun,x(n)+h/2,y(n)+h/2*k1); k3=feval(fun,x(n)+h/2,y(n)+h/2*k2); k4=feva
10、l(fun,x(n+1),y(n)+h*k3); y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6;endx=x;y=y;在MATLAB输入以下程序： clear all; fun=inline( y-x2+1); x,y= Runge_Kutta (fun,0,2,0.5,0.1); x,y plot(x,y, g+-)结果及其图象：ans = ans = 0 0.5000 0.1000 0.6574 0.2000 0.8293 0.3000 1.0151 0.4000 1.2141 0.5000 1.4256 0.6000 1.6489 0.7000 1.8831 0.8
11、000 2.1272 0.9000 2.3802 1.0000 2.6409 1.1000 2.9079 1.2000 3.1799 1.3000 3.4553 1.4000 3.7324 1.5000 4.0092 1.6000 4.2835 1.7000 4.5530 1.8000 4.8152 1.9000 5.0670 2.0000 5.30555.1.5 问题讨论与分析由以上数值分析结果绘制表格：精确解欧拉方法改进的欧拉方法四阶龙格-库塔方法xiyiyi误差yi误差yi误差00.50000.500000.500000.500000.10.65740.6500-0.00740.6570
12、-0.00040.657400.20.82930.8140-0.01530.8284-0.00090.829300.31.01510.9914-0.02371.0137-0.00141.015100.41.21411.1815-0.0326 1.2122-0.00191.214100.51.42561.3837-0.04191.4232-0.00241.425600.61.64891.5971-0.05181.6459-0.00301.648900.71.88311.8208-0.0623 1.8794-0.00371.883100.82.12722.0538-0.0734 2.1228-0.
13、00442.127200.92.38022.2952-0.0850 2.3750-0.00522.380201.02.64092.5438-0.0971 2.6348-0.00612.640901.12.90792.7981-0.10982.9010-0.00692.907901.23.17993.0569-0.12303.1720-0.00793.179901.33.45543.3186-0.13683.4464-0.00903.4553-0.00011.43.73243.5815-0.15093.7223-0.01013.732401.54.00923.8437-0.16553.9978-
14、0.01144.009201.64.28354.1030-0.18054.2708-0.01274.283501.74.55304.3573-0.19574.5390-0.01404.553001.84.81524.6040-0.21124.7996-0.01564.815201.95.06714.8405-0.22665.0499-0.01725.0670-0.00012.05.30555.0635-0.24205.2866-0.01895.30550x=0:0.1:2;y1 =0.5000 0.6574 0.8293 1.0151 1.2141 1.4256 1.6489 1.8831 2
15、.1272 2.3802 2.6409 2.9079 3.1799 3.4554 3.7324 4.0092 4.2835 4.5530 4.8152 5.0671 5.3055; y1 =0.5000 0.6574 0.8293 1.0151 1.2141 1.4256 1.6489 1.8831 2.1272 2.3802 2.6409 2.9079 3.1799 3.4554 3.7324 4.0092 4.2835 4.5530 4.8152 5.0671 5.3055; y2=0.5000 0.6500 0.8140 0.9914 1.1815 1.3837 1.5971 1.820
16、8 2.0538 2.2952 2.5438 2.7981 3.0569 3.3186 3.5815 3.8437 4.1030 4.3573 4.6040 4.840 5.0635; y3=0.5000 0.6570 0.8284 1.0137 1.2122 1.4232 1.6459 1.8794 2.1228 2.3750 2.6348 2.9010 3.1720 3.4464 3.7223 3.9978 4.2708 4.5390 4.7996 5.0499 5.2866; y4=0.5000 0.6574 0.8293 1.0151 1.2141 1.4256 1.6489 1.88
17、31 2.1272 2.3802 2.6409 2.9079 3.1799 3.4553 3.7324 4.0092 4.2835 4.5530 4.8152 5.0670 5.3055; plot(x,y1,r+-) hold on,plot(x,y2,b-) plot(x,y1,r+-) hold on,plot(x,y3,b-) plot(x,y1,r+-) hold on,plot(x,y4,b-)由上表和图可以看出欧拉法误差最大，而改进欧拉和龙格库塔方法误差相对较小，并且龙格库塔方法误差最小且大部分值都跟精确值相同。由欧拉图与精确图相比可清晰看到，随着x的增加，函数值与精确值的偏差越
18、来越大。5.2 选择用欧拉方法，改进欧拉方法，4阶龙格库塔方法之一取不同步长分别求下面微分方程的初值dy/dx=y-x2+1 y(0)=0.5 x0，2，并对结果进行分析说明，给出你的结论。用改进欧拉方法取不同步长分别求微分方程的初值dy/dx=y-x2+1 y(0)=0.5 x0，2。一步长 0.05在MATLAB输入以下程序： clear all fun=inline( y-x2+1); x,y=gjEuler(fun,0,2,0.5,0.05); x,y plot(x,y,r*-)结果及其图像：ans = 0 0.5000 0.0500 0.5768 0.1000 0.6573 0.15
19、00 0.7414 0.2000 0.8291 0.2500 0.9202 0.3000 1.0147 0.3500 1.1126 0.4000 1.2136 0.4500 1.3178 0.5000 1.4250 0.5500 1.5352 0.6000 1.6482 0.6500 1.7639 0.7000 1.8822 0.7500 2.0030 0.8000 2.1261 0.8500 2.2514 0.9000 2.3789 0.9500 2.5082 1.0000 2.6393 1.0500 2.7720 1.1000 2.9061 1.1500 3.0415 1.2000 3.1
20、779 1.2500 3.3152 1.3000 3.4531 1.3500 3.5913 1.4000 3.7298 1.4500 3.8682 1.5000 4.0063 1.5500 4.1437 1.6000 4.2803 1.6500 4.4156 1.7000 4.5494 1.7500 4.6814 1.8000 4.8112 1.8500 4.9384 1.9000 5.0627 1.9500 5.1836 2.0000 5.3007二步长0.1在MATLAB输入以下程序： clear all fun=inline( y-x2+1); x,y=gjEuler(fun,0,2,0
21、.5,0.1); x,y plot(x,y,m+-)结果及其图象： 0 0.5000 0.1000 0.6570 0.2000 0.8284 0.3000 1.0137 0.4000 1.2122 0.5000 1.4232 0.6000 1.6459 0.7000 1.8794 0.8000 2.1228 0.9000 2.3750 1.0000 2.6348 1.1000 2.9010 1.2000 3.1720 1.3000 3.4464 1.4000 3.7223 1.5000 3.9978 1.6000 4.2708 1.7000 4.5390 1.8000 4.7996 1.900
22、0 5.0499 2.0000 5.2866三步长0.15在MATLAB输入以下程序： clear all fun=inline( y-x2+1); x,y=gjEuler(fun,0,2,0.5,0.15); x,y plot(x,y,bo-)结果及其图像ans = 0 0.5000 0.1500 0.7402 0.3000 1.0121 0.4500 1.3136 0.6000 1.6422 0.7500 1.9950 0.9000 2.3687 1.0500 2.7594 1.2000 3.1625 1.3500 3.5728 1.5000 3.9842 1.6500 4.3896 1.
23、8000 4.7809 1.9500 5.1484四步长0.2在MATLAB输入以下程序： clear all fun=inline( y-x2+1); x,y=gjEuler(fun,0,2,0.5,0.2); x,y plot(x,y,yd-)结果及其图像ans = 0 0.5000 0.2000 0.8260 0.4000 1.2069 0.6000 1.6372 0.8000 2.1102 1.0000 2.6177 1.2000 3.1496 1.4000 3.6937 1.6000 4.2351 1.8000 4.7556 2.0000 5.2331五不同步长之比图像六由以上数值分
24、析结果绘制表格精确解步长为0.05步长为0.1步长为0.15步长为0.2xiyiyi误差yi误差yi误差yi误差0.0000 0.5000 0.5000 00.5000 0.0000 0.5000 0.0000 0.5000 0.0000 0.0500 0.5769 0.5768 -0.00010.1000 0.6574 0.6573 -0.00010.6570 -0.0004 0.1500 0.7416 0.7414 -0.00020.7402 -0.0014 0.2000 0.8293 0.8291 -0.00020.8284 -0.0009 0.8260 -0.0033 0.2500 0
25、.9205 0.9202 -0.00030.3000 1.0151 1.0147 -0.00041.0137 -0.0014 1.0121 -0.0030 0.3500 1.1130 1.1126 -0.00040.4000 1.2141 1.2136 -0.00051.2122 -0.0019 1.2069 -0.0072 0.4500 1.3183 1.3178 -0.00051.3136 -0.0047 0.5000 1.4256 1.4250 -0.00061.4232 -0.0024 0.5500 1.5359 1.5352 -0.00070.6000 1.6489 1.6482 -
26、0.00071.6459 -0.0030 1.6422 -0.0067 1.6372 -0.0117 0.6500 1.7647 1.7639 -0.00080.7000 1.8831 1.8822 -0.00091.8794 -0.0037 0.7500 2.0040 2.0030 -0.0011.9950 -0.0090 0.8000 2.1272 2.1261 -0.00112.1228 -0.0044 2.1102 -0.0170 0.8500 2.2527 2.2514 -0.00130.9000 2.3802 2.3789 -0.00132.3750 -0.0052 2.3687 -0.0115 0.9500 2.5096 2.5082 -0.00141.0000 2.6409 2.6393 -0.00162.6348 -0.0061 2.6177 -0.0232 1.0500 2.7737 2.7720 -0.00172.7594 -0.0143 1.1000 2.9079 2.9061 -0.00182.9010 -0.0069 1.1500 3.0434 3.0415
