《常微分方程数值解法课程设计报告1.doc》由会员分享，可在线阅读，更多相关《常微分方程数值解法课程设计报告1.doc（14页珍藏版）》请在三一文库上搜索。
1、数值分析课程设计数值分析课程设计题目：常微分方程的数值解法题目：常微分方程的数值解法 组员：邓全飞组员：邓全飞 201120320105 许曦许曦 201120320108 熊鲁平熊鲁平 201120320203 余佳明余佳明 201120320206 专业：信息与计算科学专业：信息与计算科学 指导老师：刘唐伟指导老师：刘唐伟 一、摘要一、摘要.3二、设计目的二、设计目的 .3三、理论基础三、理论基础 .31.1.欧拉公式欧拉公式【3】【3】.32.2.改进改进 EulerEuler 方法方法【3】【3】.33.3.三阶龙格三阶龙格- -库塔方法库塔方法【3】【3】.34.4.四阶龙格四阶龙格
2、- -库塔方法库塔方法【3】【3】.4四、程序代码及运算结果四、程序代码及运算结果 .41 1、用欧拉法求解、用欧拉法求解 .52 2、用改进欧拉法求解：、用改进欧拉法求解： .63 3、用、用 3 3 阶龙格阶龙格库塔求解库塔求解 .74 4、 用用 4 4 阶龙格阶龙格库塔求解库塔求解 .96 6、欧拉方法与改进欧拉方法、欧拉方法与改进欧拉方法、3 3 阶龙格库塔法以及阶龙格库塔法以及 4 4 阶龙格库塔法得出的解得比较。阶龙格库塔法得出的解得比较。 .10五、数值分析设计的五、数值分析设计的 GUIGUI 界面界面.13六、结果分析六、结果分析 .14七、设计心得七、设计心得 .14八、
3、参考文献八、参考文献 .14一、摘要一、摘要在 matlab 环境下熟悉的运用计算机编程语言并结合龙格库塔法的理论基础对常微分方程初值问题进行求解，在运行完程序后以及对运行结果做出各方面的分析和比较。二、设计目的二、设计目的用熟悉的计算机语言编程上机完成用欧拉方法、改进欧拉方法、3 阶龙格库塔法以及 4 阶龙格库塔法求解常微分方程初值问题。三、理论基础三、理论基础1.1.欧拉公式欧拉公式【3】【3】在点 将作 Taylor 展开，得nx1()()nny xy xh，那么当 h 充分小时，略211()()()(),(,)2!nnnnnnnhy xy xhy xyxx去误差项，用近似替代、近似替代
4、，并注意到2()2!nhyny()ny x1ny1()ny x，便得()(, ()nnny xf xy x0010(),(,),0,1,1,.nnnnnyy xyyhf xyxxnhbanNhN 上述方法称为 Euler 方法。2.2.改进改进 EulerEuler 方法方法【3】【3】在应用梯形方法的迭代公式进行运算时，每迭代一次都要重新计算函数的值，且还要判断何时可以终止或转下一步计算。为了控制计算量和简( , )f x y化算法，通常只迭代一次就转入下一步计算。具体说，我们先用 Euler 公式求得一个初步的近似值，称之为预测值，然后用梯形方法的迭代公式作一次1ny迭代得，即将校正一次，
5、这样建立的预测-校正方法称之为改进的1ny1nyEuler 方法：预测：，1(,)nnnnyyhf xy校正：. .111 (,)(,)2nnnnnnhyyf xyf xy3.3.三阶龙格三阶龙格- -库塔方法库塔方法【3】【3】类似前面改进的 Euler 方法公式的推导方法，将在处作 Taylor1ny(,)nnxy展开，然后再将在处作 Taylor 展开，只要将两个展开式前四项相1()ny xnxx同便有。于是得到三阶龙格-库塔公式为：411()()nny xyo h4 4.四阶龙格四阶龙格- -库塔方法库塔方法【3】【3】 类似前面三阶龙格-库塔的推导方法，如果每步计算四次函数 f（x，
6、y）的值，完全类似的，可以导出局部截断误差为的四阶龙格-库塔5()o h公式，其公式为：112341213243(22),6(,),(,),2211(,),22(,).nnnnnnnnnnhyykkkkkf xyhhkf xykkf xh yhkkf xh yhk四、程序代码及运算结果四、程序代码及运算结果用欧拉方法、改进欧拉方法、3 阶龙格库塔法以及 4 阶龙格库塔法求解常微分方程初值问题：2(0)2dyxyydx 精确解为：y0= -2.0000 -1.9802 -1.9231 -1.8349 -1.7241 -1.6000 -1.4706 -1.3423 -1.2195 -1.1050
7、-1.0000 1123121312(4),6(,),(,),22(,2).nnnnnnnnhyykkkkf xyhhkf xykkf xh yhkhk 精确解图形1、用欧拉法求解程序如下：建立函数文件 cwfa1.mfunction x,y=cwfa1(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( y2*x ); x,y=cwfa1(fu
8、n,0,1,-2,0.1); x,y plot(x,y,r+-)ans = 0 -2.0000 0.1000 -2.0000 0.2000 -1.9600 0.3000 -1.8832 0.4000 -1.7768 0.5000 -1.6505 0.6000 -1.5143 0.7000 -1.3767 0.8000 -1.2440 0.9000 -1.1202 1.0000 -1.0073结果及其图象： 图 1 欧拉法解 2、用改进欧拉法求解：程序如下：建立函数文件 cwfa2.mfunction x,y=cwfa2(fun,x_span,y0,h)x=x_span(1):h:x_span(
9、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 fun=inline( y2*x); x,y=cwfa2(fun,0,1,-2,0.1); x,y plot(x,y,r+-) ans = 0 -2.0000 0.1000 -1.9800 0.2000 -1.9227 0.3000 -1.8345 0.4000 -1.72
10、39 0.5000 -1.6001 0.6000 -1.4711 0.7000 -1.3432 0.8000 -1.2208 0.9000 -1.10661.0000 -1.0018结果及其图像： 图 3 改进欧拉解3、用 3 阶龙格库塔求解程序如下：建立函数文件 cwfa4.mfunction x,y=cwfa4(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,
11、x(n)+h,y(n)+h*(2*k2-k1); y(n+1)=y(n)+h*(k1+4*k2+k3)/6;endx=x;y=y;在 Matlab 输入以下程序：clear all;fun=inline( y2*x);x,y=cwfa4(fun,0,1,-2 ,0.1);x,y plot(x,y, b*-)ans = 0 -2.0000 0.1000 -1.9803 0.2000 -1.9232 0.3000 -1.8350 0.4000 -1.7243 0.5000 -1.6002 0.6000 -1.4707 0.7000 -1.3424 0.8000 -1.2196 0.9000 -1.
12、1050 1.0000 -1.0000结果及图像： 图 3 3 阶龙格-库塔解4、 用 4 阶龙格库塔求解程序如下：建立函数文件 cwfa3.mfunction x,y=cwfa3(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=feval(fun,x(n+1),y(n)+h*k3); y(n+1)=y(n)+h
13、*(k1+2*k2+2*k3+k4)/6;endx=x;y=y;在 MATLAB 输入以下程序： clear all; fun=inline( y2*x ); x,y=cwfa3(fun,0,1,-2 ,0.1); x,y plot(x,y, b*-) ans = 0 -2.0000 0.1000 -1.9802 0.2000 -1.9231 0.3000 -1.8349 0.4000 -1.7241 0.5000 -1.6000 0.6000 -1.4706 0.7000 -1.3423 0.8000 -1.2195 0.9000 -1.1050 1.0000 -1.0000结果及其图象：
14、图 4 4 阶龙格-库塔解6、欧拉方法与改进欧拉方法、3 阶龙格库塔法以及 4 阶龙格库塔法得出的解得比较。x=0:0.1:1;y0= -2.0000 -1.9802 -1.9231 -1.8349 -1.7241 -1.6000 -1.4706 -1.3423 -1.2195 -1.1050 -1.0000y1= -2.0000 -2.0000 -1.9600 -1.8832 -1.7768 -1.6505 -1.5143 -1.3767 -1.2440 -1.1202 -1.0073;y2= -2.0000 -1.9800 -1.9227 -1.8345 -1.7239 -1.6001 -
15、1.4711 -1.3432 -1.2208 -1.1066 -1.0018;y3= -2.0000 -1.9803 -1.9232 -1.8350 -1.7243 -1.6002 -1.4707 -1.3424 -1.2196 -1.1050 -1.0000;y4= -2.0000 -1.9802 -1.9231 -1.8349 -1.7241 -1.6000 -1.4706 -1.3423 -1.2195 -1.1005 -1.1000;plot(x,y0,r*-)hold on,plot(x,y1,b-),title(精确解与欧拉方法比较)plot(x,y0,r+-)hold on,pl
16、ot(x,y2,b-),title(精确解与改进欧拉方法比较)plot(x,y0,r+-)hold on,plot(x,y3,b-),title(精确解与 3 阶龙格库塔方法比较)plot(x,y0,r+-)hold on,plot(x,y4,b-),title(精确解与 4 阶龙格库塔方法比较) 图 5 精确解与欧拉法比较 图 6 精确解与改进欧拉法比较 图 7 精确解与 3 阶-龙格库塔方法比较图 8 精确解与 4 阶龙格-库塔方法比较五、数值分析设计的五、数值分析设计的 guigui 界面界面六、结果分析六、结果分析由以上的各种方法与精确解的比较图可以看出在计算精度上，四阶龙格-库塔方法
17、的误差最小，其次是三阶龙格-库塔方法，再次之就是改进欧拉方法，欧拉方法误差则比较大，所以四阶龙格-库塔方法得到最佳的精度。而在计算量上面，相应地，很明显的四阶龙格-库塔方法也是最大，三阶龙格-库塔次之，其后为改进欧拉方法，欧拉方法计算量最小。这样的结果，说明了运用以上三种方法时，其计算量的多少与精度的大小成正比。我们在实际运用与操作中，可以根据实际情况，选择这 4 种方法中的其中一种最适合的，追求精度的话，可以使用三阶或四阶龙格-库塔方法；而改进的欧拉方法，在精度上和计算量上都表现得很出色，能够满足一般情况；而欧拉方法更主要的是适用于对的估计ny上，相应的，精度则有所欠缺。以上的选择，都取决于
18、具体的情况。七、设计心得七、设计心得 这次实验花了较多的时间，先是选择方程方面遇到了困难：有些方程在运行过程中会不断出错，错误的原因有很多，有的是在输入格式上有错，有的是在步长和初值的取法上出错。在试了十多个方程之后终于找到一个自己满意的。接下来的运行中也出现了一些错误。在用第一个方程求精确解时由于没有区分点乘和乘导致错误，有些方程精确解时可以运行但是在用欧拉方程和改进的欧拉方法时就会报错，有的方程即使求的出来结果与精确解的误差是很大的，显然是错误的。还有就是刚开始犯了一个低级错误：把原函数直接代入欧拉求解这样当然是不正确的。取初值对函数图像的影响很大，如果初值取的不合适误差也会非常大，所以要合理的选择初值，在几种求解方法中只有龙格库塔的求解与精确值最为接近几乎是重合的。通过的后面的不同步长求解的比较发现步长越小与精确值越接近。总之在做这次实验中对学习数学有了更多的体会，那就是一定要认真，不能丝毫马虎，一个小小的错误就会导致所有的数据不能运行。八、参考文献八、参考文献1 刘卫国.MATLAB 程序设计与应用程. 北京: 高等教育出版社，20092奚如成 数值分析方法 中国科学技术大学出版社。