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

常微分方程的数值解法 专业班级：信息软件 姓名：吴中原 学号：120108010002 一、实验目的 2、明确各种算法的精度与所选步长有密切关系；通过计算更加了解各种 算法的优越性。 二、实验题目 2、试分别取不同步长，考察某节点处数值解的误差变化情况； 3、试用不同算法求解某初值问题，结果有何异常； 4、分析各个算法的优缺点。 三、实验原理与理论基础 (6-1) (6-2) 用数值方法求解时，我们总是认为(6-1)、(6-2)的解存在且唯一。 欧拉法是解初值问题的最简单的数值方法。从(6-2)式由于y (x0) = y0已给定，因而可以算出 设x1 = h充分小，则近似地有： (6-3) 记 从而我们可以取 作为的近似值。利用及f (x1, y1)又可以算出的近似值： 一般地，在任意点处的近似值由下式给出 (6-4) 这就是欧拉法的计算公式，h称为步长。 ，它每一步计算的值一次，截断误差为。 改进的欧拉公式可以改写为： ，它每一步要计算的值两次，截断误差为。 改进的欧拉方法之所以比欧拉方法具有更高的精度，是因为在每一步它都比欧拉方法多计算了一次的值。因此，要进一步提高精度，可以考虑在每一步增加计算的次数。 如果考虑在每一步计算的值四次，则可以推得如下公式： 此公式称为标准四阶龙格-库塔(Runge-Kutta)公式，它的截断误差为。虽然用龙格-库塔方法每一步需要四次调用，计算量较改进的欧拉方法大一倍，这里由于龙格-库塔方法的步长增大了一倍，因而两种方法总的计算量相同，但龙格-库塔方法精确度更高。所以龙格-库塔公式兼顾了精度和计算工作量的较为理想的公式，在实际计算中最为常用。 实验内容 科学计算中经常遇到微分方程（组）初值问题，需要利用Euler法，改进Euler 法，Rung-Kutta方法求其数值解，诸如以下问题： （1）0 x1 分别取h=0.1,0.2,0.4时数值解。 初值问题的精确解。 （2）r=3的Adams显式和预 - 校式求解 取步长h=0.1，用四阶标准R-K方法求值。 （3）Euler法或四阶标准R-K方法求解 取步长0.01,计算数值解，参考结果 （4）R- K方法求二阶方程初值问题的数值解 （I） (II) (III) (IV) 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=1:M y(k+1)=y(k)+h*t(k)*y(k)^(1/3); end yb=y(M+1); yy=((t.^2+2)./3).^1.5; det=yy-y; plot(t,y,r-,t,yy,b:,t,det); 改进欧拉法程序 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+1),y(k)+h*p); y(k+1)=y(k)+0.5*h*(p+q); end yy=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;K2=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);