常微分方程的数值解法实验报告

1、精选优质文档-倾情为你奉上常微分方程的数值解法专业班级：信息软件 姓名：吴中原 学号：2 一、实验目的 1、熟悉各种初值问题的算法，编出算法程序； 2、明确各种算法的精度与所选步长有密切关系；通过计算更加了解各种 算法的优越性。 二、实验题目 1、根据初值问题数值算法，分别选择二个初值问题编程计算； 2、试分别取不同步长，考察某节点处数值解的误差变化情况； 3、试用不同算法求解某初值问题，结果有何异常； 4、分析各个算法的优缺点。三、实验原理与理论基础（一） 欧拉法算法设计对常微分方程初始问题 (6-1) (6-2)用数值方法求解时，我们总是认为(6-1)、(6-2)的解存在且唯一。欧拉法是解

2、初值问题的最简单的数值方法。从(6-2)式由于y (x0) = y0已给定，因而可以算出。设x1 = h充分小，则近似地有：(6-3)记 从而我们可以取作为的近似值。利用及f (x1, y1)又可以算出的近似值：一般地，在任意点处的近似值由下式给出(6-4)这就是欧拉法的计算公式，h称为步长。（2） 四阶龙格-库塔法算法设计：欧拉公式可以改写为：，它每一步计算的值一次，截断误差为。改进的欧拉公式可以改写为： ，它每一步要计算的值两次，截断误差为。改进的欧拉方法之所以比欧拉方法具有更高的精度，是因为在每一步它都比欧拉方法多计算了一次的值。因此，要进一步提高精度，可以考虑在每一步增加计算的次数。如

3、果考虑在每一步计算的值四次，则可以推得如下公式：此公式称为标准四阶龙格-库塔(Runge-Kutta)公式，它的截断误差为。虽然用龙格-库塔方法每一步需要四次调用，计算量较改进的欧拉方法大一倍，这里由于龙格-库塔方法的步长增大了一倍，因而两种方法总的计算量相同，但龙格-库塔方法精确度更高。所以龙格-库塔公式兼顾了精度和计算工作量的较为理想的公式，在实际计算中最为常用。四、实验内容（一）问题重述： 科学计算中经常遇到微分方程（组）初值问题，需要利用Euler法，改进Euler 法，Rung-Kutta方法求其数值解，诸如以下问题：（1） 0 < x1 分别取h=0.1,0.2,0.4时数值

4、解。 初值问题的精确解。（2）用r=3的Adams显式和预 - 校式求解 取步长h=0.1，用四阶标准R-K方法求值。（3）用改进Euler法或四阶标准R-K方法求解 取步长0.01,计算数值解，参考结果 （4）利用四阶标准R- K方法求二阶方程初值问题的数值解 （I） (II) (III) (IV) (二)实验代码：1、欧拉法程序function y=Euler(a,b,M,y0)%a=1,b=2,M=10,f=t*y(1/3),y0=1;h=(b-a)/M;t=zeros(1,M+1);t=a:h:b;y=zeros(1,M+1);yy=zeros(1,M+1);y(1)=y0;for k

5、=1:M y(k+1)=y(k)+h*t(k)*y(k)(1/3);endyb=y(M+1);yy=(t.2+2)./3).1.5;det=yy-y;plot(t,y,'r-',t,yy,'b:',t,det);2、 改进欧拉法程序function H=heeuler(a,b,M,ya,f)%a=0,b=1,M=10,f=t*t+t-y,y0=0;h=(b-a)/M;t=zeros(1,M+1);y=zeros(1,M+1);p=0;q=0;t=a:h:b;y(1)=ya;for k=1:M p=feval(f,t(k),y(k); q=feval(f,t(k+

6、1),y(k)+h*p); y(k+1)=y(k)+0.5*h*(p+q);endyy=t.*t-t+1-exp(-t);det=yy-y;plot(t,y,'r-',t,yy,'b:',t,det);H=t',y',yy',det'function f=ff(t,y);f=t.2+t-y;3、四阶龙格-库塔法程序function H=r_k4(a,b,M,ya,f)%a=0,b=1,M=10,f=t*t+t-y,y0=0;h=(b-a)/M;t=zeros(1,M+1);t=a:h:b;y=zeros(1,M+1);K1=0;K

7、2=0;K3=0;K4=0;y(1)=ya;for k=1:M K1=feval(f,t(k),y(k); K2=feval(f,t(k)+0.5*h,y(k)+0.5*h*K1); K3=feval(f,t(k)+0.5*h,y(k)+0.5*h*K2); K4=feval(f,t(k)+h,y(k)+h*K3); y(k+1)=y(k)+1/6*(K1+2*K2+2*K3+K4);endyy=t.*t-t+1-exp(-t);det=yy-y;plot(t,y,t,yy,t,det);H=t',y',yy',det'function f=ff(t,y);f=

8、t.2+t-y; 五、实验结果1） 0 < x1 分别取h=0.1,0.2,0.4时数值解。 初值问题的精确解。Euler('han',0,0,2,0.1)ans = 0 0 0.1000 0 0.2000 0.0400 0.3000 0.1192 0.4000 0.2356 0.5000 0.3862 0.6000 0.5669 0.7000 0.7729 0.8000 0.9988 0.9000 1.2389 1.0000 1.4874 1.1000 1.7386 1.2000 1.9874 1.3000 2.2289 1.4000 2.4591 1.5000 2.6

9、749 1.6000 2.8736 1.7000 3.0539 1.8000 3.2147 1.9000 3.3561 2.0000 3.4784Euler('han',0,0,2,0.2)ans = 0 0 0.2000 0 0.4000 0.1600 0.6000 0.4672 0.8000 0.8911 1.0000 1.3886 1.2000 1.9108 1.4000 2.4122 1.6000 2.8568 1.8000 3.2226 2.0000 3.5025Euler('han',0,0,2,0.4)ans = 0 0 0.4000 0 0.8000 0.6400 1.2000 1.7152 1.6000 2.8119 2.0000 3.5723 2、四阶龙格-库塔法结果 取步长h=0.1，用四阶标准R-K方法求值。RK4('han',-1,0,1,0.1)ans = -1.0000 1.0000 -0.9000 0.9909 -0.8000 0.9672 -0.7000 0.9331 -0.6000 0.8921 -0.5000 0.8468 -0.4000 0.7993 -0.3000 0.7515 -0.2000 0.7048 -0.1000