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

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

2、法是解初值问题的最简单的数值方法。从(6-2)式由于y(x0)=y0已给定,因而可以算出y'(x)二f(x,y)。000设x1=h充分小，则近似地有：y0)(6-3),ni记y=y(x)i=0,1,ii从而我们可以取y二y二hf(x,y)1000作为y(xi)的近似值。利用y1及fg,yj又可以算出y(x2)的近似值:y2=yi+hf(片,yi)一般地，在任意点x=(n+1)h处y(x)的近似值由下式给出n+1(6-4)y=y+hf(x,y)n+1nnn这就是欧拉法的计算公式，h称为步长。二)四阶龙格-库塔法算法设计：欧拉公式可以改写为：广1甘(1)，匕每步计算f(x,y)的值次，截1

3、1ii断误差为o)。y=y+(k+k)i+1i22改进的欧拉公式可以改写为：<k=hf(x,y)，它每一步要计算k=hf(x+h,y+k)2ii1f(x,y)的值两次，截断误差为o(h3)。改进的欧拉方法之所以比欧拉方法具有更高的精度，是因为在每一步它都比欧拉方法多计算了一次f(x,y)的值。因此，要进一步提高精度，可以考虑在每一步增加计算f(x,y)的次数。如果考虑在每一步计算f(x,y)的值四次，则可以推得如下公式：y=y+1(k+2k+2k(xi+1i61231iih1k=hfx+,y+k2Ij2i21丿(h1k=hfx+,y+k3Ji2i22丿k=hf(x,y)k4hfi+h,y

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

5、数值解。初值问题的精确解y7+5e-x2。用r=3的Adams显式和预-校式求解y'=X2-y23)、y(-1)=01<x<0取步长h=0.1,用四阶标准R-K方法求值。用改进Euler法或四阶标准R-K方法求解y1G)=-1y2G)=0y3G)=10<x<1(、0.01,计算y血)，血10)；血15)数值解，参y615)q0.9880787,y(0.15)q0.1493359,y(0.15人0.8613125123(4)利用四阶标准R-K方法求二阶方程初值问题的数值解3y，+2y=0=0,y，G)=10<x<1h=0.02y1=y2<y2=y

6、1、y3=-y3取步长h=0.9880787,y考结果I)(II)y-0.1(1-y2)y'+y=0、yG)=1,y(0)=00<x<1,h=0.1(III)y'二-7<ex+1y(0)=1,y心)=00<x<2,h=0.1y"+siny=0(IV)M)=i,y心)=0(二)实验代码：1、欧拉法程序functiony=Euler(a,b,M,y0)%a=1,b=2,M=10,f=t*yA(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(

7、1)=y0;fork=1:My(k+1)=y(k)+h*t(k)*y(k)人(1/3);endyb=y(M+1);0<x<4,h=0.2yy=(t.A2+2)./3)41.5;det=yy-y;plot(t,y'r-',t,yy,'b:',t,det);2、改进欧拉法程序functionH=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;fork=1:Mp=feval(f,t(

8、k),y(k);q=feval(f,t(k+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'functionf=ff(t,y);f=t.入2+t-y;3、四阶龙格-库塔法程序functionH=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=

9、zeros(1,M+1);K1=0;K2=0;K3=0;K4=0;y(1)=ya;fork=1:MK1=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

10、f=ff(t,y);f=t.入2+t-y;五、实验结果,4xy=-xy<y1)yG)=0分别取h=0.1,0.2,0.4时数值解。初值问题的精确解y='''4+5e-x2oEuler('han',0,0,2,0.1)ans=000.100000.20000.04000.30000.11920.40000.23560.50000.38620.60000.56690.70000.77290.80000.99880.90001.23891.00001.48741.10001.73861.20001.98741.30002.22891.40002.4591

11、1.50002.67491.60002.87361.70003.05391.80003.21471.90003.35612.00003.4784Euler('han',0,0,2,0.2)ans=000.200000.40000.16000.60000.46720.80000.89111.00001.38861.20001.91081.40002.41221.60002.85681.80003.22262.00003.5025Euler('han',0,0,2,0.4)ans=000.400000.80000.64001.20001.71521.60002.81192.00003.57232、四阶龙格-库塔法结果y'=X2-y2V、y(1L0-1<x<0取步长h=0.1,用四阶标准R-K方法求值。RK4('han',-1,0,1,0.1)ans=-1.00001.0000-0.90000.9909-0.80000.9672-0.70000.9331-0.60000.8921-0.50000.8468-0.