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

常微分方程的数值解法专业班级：信息软件姓名：吴中原学号：120108010002一、实验目的1、熟悉各种初值问题的算法，编出算法程序；2、明确各种算法的精度与所选步长有密切关系；通过计算更加了解各种算法的优越性。二、实验题目1、根据初值问题数值算法，分别选择二个初值问题编程计算；2、试分别取不同步长，考察某节点jx处数值解的误差变化情况；3、试用不同算法求解某初值问题，结果有何异常；4、分析各个算法的优缺点。三、实验原理与理论基础（一）欧拉法算法设计对常微分方程初始问题(6-1)(6-2)用数值方法求解时，我们总是认为(6-1)、(6-2)的解存在且唯一。欧拉法是解初值问题的最简单的数值方法。从(6-2)式由于y(x0)=y0已给定，因而可以算出),()('000yxfxy。设x1=h充分小，则近似地有：),()(')()(00001yxfxyhxyxy(6-3)记,n,,ixyyii10)(从而我们可以取),(0001yxhfyy作为)(1xy的近似值。利用1y及f(x1,y1)又可以算出)(2xy的近似值：),(1112yxhfyy一般地，在任意点hnxn11处)(xy的近似值由下式给出),(1nnnnyxhfyy(6-4)这就是欧拉法的计算公式，h称为步长。)(),(dd00yxyyxfxy（二）四阶龙格-库塔法算法设计：欧拉公式可以改写为：111,iiiiyykkhfxy，它每一步计算,fxy的值一次，截断误差为2oh。改进的欧拉公式可以改写为：11212112,,iiiiiiyykkkhfxykhfxhyk，它每一步要计算,fxy的值两次，截断误差为3oh。改进的欧拉方法之所以比欧拉方法具有更高的精度，是因为在每一步它都比欧拉方法多计算了一次,fxy的值。因此，要进一步提高精度，可以考虑在每一步增加计算,fxy的次数。如果考虑在每一步计算,fxy的值四次，则可以推得如下公式：1123412132431226,1,221,22,iiiiiiiiiiyykkkkkhfxyhkhfxykhkhfxykkhfxhyk此公式称为标准四阶龙格-库塔(Runge-Kutta)公式，它的截断误差为5oh。虽然用龙格-库塔方法每一步需要四次调用f，计算量较改进的欧拉方法大一倍，这里由于龙格-库塔方法的步长增大了一倍，因而两种方法总的计算量相同，但龙格-库塔方法精确度更高。所以龙格-库塔公式兼顾了精度和计算工作量的较为理想的公式，在实际计算中最为常用。四、实验内容（一）问题重述：科学计算中经常遇到微分方程（组）初值问题，需要利用Euler法，改进Euler法，Rung-Kutta方法求其数值解，诸如以下问题：（1）0yxyyxy040x1分别取h=0.1,0.2,0.4时数值解。初值问题的精确解254xey。（2）用r=3的Adams显式和预-校式求解2201yxyy01x取步长h=0.1，用四阶标准R-K方法求值。（3）用改进Euler法或四阶标准R-K方法求解10330012101213yyy2yyyyyy10x取步长h0.01,计算15.0,10.0,05.0yyy数值解，参考结果8613125.015.0,1493359.015.0,9880787.015.0321yyy（4）利用四阶标准R-K方法求二阶方程初值问题的数值解（I）10,00023yyyyy0.02hx10(II)00,100)1(1.02yyyyyy0.1hx,10(III)00,101yyeyyx0.1hx,20(IV)00,100sinyyyy2,400.hx(二)实验代码：1、欧拉法程序functiony=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;fork=1:My(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、改进欧拉法程序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(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=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']functionf=ff(t,y);f=t.^2+t-y;五、实验结果1）0yxyyxy040x1分别取h=0.1,0.2,0.4时数值解。初值问题的精确解254xey。Euler('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.45911.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、四阶龙格-库塔法结果2201yxyy01x取步长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.40000.7993-0.30000.7515-0.20000.7048-0.10000.660600.6199六、实验结果分析与小结1、步长h越小则计算精度越高，相对误差越小。因此，在计算能力允许的范围内，步长越小，得到的结果更精确。2、改进欧拉法和欧拉法相比较，改进欧拉法的计算精度要更高，相对误差也较小。因此在求常微分方程的数值解时，改进欧拉法比欧拉法更精确。