(2)常微分方程的差分方法.doc

实验三：常微分方程的差分方法程序设计与验证李郭 应数081 0809501024一、 实验结果1. 改进的欧拉方法（1）建立euler函数文件(euler.m)euler.m文件程序如下：functionx,y=euler(dyfun,xspan,y0,h)% dyfun为函数f(x,y)%xspan为求解区间x0,xN%h为步长%x返回节点%y返回数值解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）在工作窗口输入如下程序ticformat longdyfun=inline(y-2*x/y);x,y=euler(dyfun,0,1,1,0.2);x,ytoc（3）结果步长时，结果为： 0 1.00000.1000 1.09590.2000 1.18410.3000 1.26620.4000 1.34340.5000 1.41640.6000 1.48600.7000 1.55250.8000 1.61650.9000 1.67821.0000 1.7379步长时，结果为： 0 1.0000 0.2000 1.1867 0.4000 1.3483 0.6000 1.4937 0.8000 1.6279 1.0000 1.75422. 四阶经典龙格-库塔方法(1)建立nark4函数文件nark4.m程序如下functionx,y=nark4(dyfun,xspan,y0,h)% dyfun为函数f(x,y)%xspan为求解区间x0,xN%h为步长%x返回节点%y返回数值解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,x(n)+h/2,y(n)+h/2*k1); k3=feval(dyfun,x(n)+h/2,y(n)+h/2*k2); k4=feval(dyfun,x(n+1),y(n)+h*k3); y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6;endx=x；y=y;（2）在工作窗口输入如下程序dyfun=inline(y-2*x/y);x,y=nark4(dyfun,0,1,1,0.1);x,y（3）结果步长时，结果为： 0 1.00000.1000 1.09540.2000 1.18320.3000 1.26490.4000 1.34160.5000 1.41420.6000 1.48320.7000 1.54920.8000 1.61250.9000 1.67331.0000 1.7321步长时，结果为： 0 1.00000.2000 1.18320.4000 1.34170.6000 1.48330.8000 1.61251.0000 1.7321二、 分析讨论 由结果可知，改进的欧拉公式和四阶经典龙格库塔方法均能求解常微分方程的初值问题，对于改进的欧拉格式，得到的结果比更精确，同时，四阶经典龙格库塔方法也是如此，且对于实验所给的题目，时，四阶经典龙格库塔方法得到的结果和真实值最为接近。aa=1.73205080756888;plot(x1,y1,r,x2,y2,b)hold onplot(aa,mh)legend(euler,nark4,精确值)grid on图一：时，两种差分方法结果与真