计算方法--常微分方程求解实验.doc

实验五 常微分方程求解实验一、 实验目的通过本实验学会对给定初值我呢他，用欧拉法、改进欧拉法、四阶龙格-库塔法求数值解和误差，并比较优缺点.对给定刚性微分方程，求其数值解，并与精确解比较，分析计算结果.二、 实验题目1. 解初值问题各种方法比较实验题目：给定初值问题 精确解为，按（1） 欧拉法，步长（2） 改进欧拉法，步长（3） 四阶标准龙格-库塔法，步长求在节点处的数值解及误差，比较各方法的优缺点.2. 刚性方程计算实验题目：给定刚性微分方程其精确解为.任选取一显示方法，取不同的步长求解，并分析计算结果.三、 实验原理1.欧拉格式由数值微分的向前差商公式可以解决初值问题（6.1）中的导数的数值计算问题：由此可得（6.1）实际上给出于是有再由得 （6.2）递推公式（6.2）称为欧拉格式。2.改进欧拉格式先对欧拉格式（6.2）对进行计算，并将结果记为，再代入（6.7）可得“预报-校正”形式的差分格式：公式（6.8）称为改进欧拉格式。3.四阶经典龙格-库塔格式：四、 实验内容1. 解初值问题各种方法比较实验程序：function shiyan51h=0.1;dyfun=inline(y./x+x*exp(x);x,y1=maeuler1(dyfun,1,2,0,h);x,y2=maeuler(dyfun,1,2,0,h);x,y3=marunge4(dyfun,1,2,0,h);y=x.*(exp(x)-exp(1);err1=abs(y-y1);err2=abs(y-y2);err3=abs(y-y3);x,err1,err2,err33.刚性方程计算实验程序：function shiyan53clc%format longfun=inline(-600*y+1199.8*exp(-0.1*x)-600);tic;x,y1=marunge4(fun,0,5,2,0.001);yb1=exp(-600*x)+2*exp(-0.1*x)-1;t1=tocerr1=abs(yb1-y1);s1=sum(err1.2)%tic;x,y2=marunge4(fun,0,5,2,0.0001);Yb2=exp(-600*x)+2*exp(-0.1*x)-1;t2=tocerr2=abs(yb2-y2);s2=sum(err2.2)五、 实验结果1. 解初值问题各种方法比较的实验结果1.00000000000000 0 0 0 1.10000000000000 0.04264443219022 0.00097356580312 0.00000092300897 1.20000000000000 0.09520401466700 0.00199790937570 0.00000178250161 1.30000000000000 0.15865732017504 0.00306315484526 0.00000259633614 1.40000000000000 0.23411778154437 0.00415909973975 0.00000337650650 1.50000000000000 0.32284614010794 0.00527488906625 0.00000413124347 1.60000000000000 0.42626522205571 0.00639872308444 0.00000486627011 1.70000000000000 0.54597705499551 0.00751757761384 0.00000558557848 1.80000000000000 0.68378242658891 0.00861692278559 0.00000629192402 1.90000000000000 0.84170305352887 0.00968043001901 0.00000698714712 2.00000000000000 1.02200658351062 0.01068965912814 0.000007672385813.刚性方程计算的实验结果t1 =8.9953s1 =1.3221e-006t2 =100.9231s2 =5.3715e-014六、 实验结果分析从实验第一小题的实验结果可以看出，用欧拉法对初值问题的求解误差相对比较大，改进欧拉法对初值问题的求解次之，而四阶龙格-库塔法对初值问题的求解误差最小.因此对初值问题的求解用四阶龙格-库塔法来求解比较准确.从实验第二小题大实验结果可以看出，当取步长比较小时，程序运行时间反而增大了，对应当误差也变大了很多，因此对给定刚性微分方程的求解，步长不应取过小.附：1.欧拉法的实验程序：%maeuler.mfunction x, y=maeuler(dyfun,xspan,y0,h)% y=f(x,y), y(x0)=y0% x,y=maeuler(dyfun,xspan,y0,h) dyfunf(x,y), xspan% x0,xn, y0y(x0), h, xy%format short;x=xspan(1):h:xspan(2); y(1)=y0;for n=1:(length(x)-1) k1=feval(dyfun,x(n),y(n); y(n+1)=y(n)+h*k1; k2=feval(dyfun,x(n+1),y(n+1); y(n+1)=y(n)+h*(k1+k2)/2;endx=x; y=y;2.改进欧拉法的实验程序：%maeuler1.mfunction x, y=maeuler1(dyfun,xspan,y0,h)% y=f(x,y), y(x0)=y0% x,y=maeuler(dyfun,xspan,y0,h) dyfunf(x,y), xspan% x0,xn, y0y(x0), h, xy%format short;x=xspan(1):h:xspan(2); y(1)=y0;for n=1:(length(x)-1) k1=feval(dyfun,x(n),y(n); y(n+1)=y(n)+h*k1;endx=x; y=y;3.四阶龙格-库塔法的实验程序：%marunge4.mfunction x,y=marunge4(dyfun,xspan,y0,h)%4y=f(x, y), y(x0)=y0%x, y=marunge4(dyfun,xspan,y0,h) dyfunf(x,y), %xspanx0, xn, y0, h, x, y %format short; x=xspan(1):h:xspan(2); y(1)=y0; for n=1:(length(x)-1) k1=feval(dyfun, x(n), y(n); k2=feval(dyfun,