资源描述：
《(2)常微分方程的差分方法》由会员上传分享，免费在线阅读，更多相关内容在工程资料-天天文库。
1、实验三：常微分方程的差分方法程序设计与验证0809501024李郭应数081一、实验结果1.改进的欧拉方法(1)建立culcr函数文件(culcr.m)euler.m文件程序如下：function[x,y]=euler(dyfun,xspan,yO,h)%dyfun为函数f(x,y)%xspan为求解区间[xO,xN]%h为步长%x返回节点%y返回数值解x=xspan(l):h:xspan(2);y(i)=y0；forn=l:length(x)-lk1=feval(dyfun,x(n),y(n));y(n+l)=y(n)+h*kl
2、;k2=feval(dyfun,x(n+1),y(n+1));y(n+1)=y(n)+h*(k1+k2)/2;endx=x';y=y：(2)在工作窗口输入如下程序ticformatlongdyfun=inline(,y-2*x/y,);[x,y]=euler(dyfun,[0,1],1,0.2);[x,y]toe(3)结果步长h=0」时，结果为:1.00001.09591.18411.266200.10000.20000.30000.40000.50000.60000.70000.80000.90001.00001.34341.4
3、1641.48601.55251.61651.67821.7379步长力=0・2时，结果为:00.20000.40000.60000.80001.00001.00001.18671.34831.49371.62791.75421.四阶经典龙格-库塔方法(1)建立nark4函数文件nark4.m程序如卜function[x,y]=nark4(dyfun,xspan,y0,h)%dyfun为函数f(x,y)%xspan为求解区间fxO,xN]%h为步长%x返回节点%y返回数值解x=xspan(l):h:xspan(2);y(i)=y0
4、；forn=l:length(x)-lkl=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+l)=y(n)+h*(kl+2*k2+2*k3+k4)/6;endx=x'；y=y：(2)在工作窗口输入如下程序dyfun=inline(,y-2*x/y,);[x,刃=nark4(dyfun,[0,1]J,0.1);[x,
5、y](1)结果步长h=0」时，结果为:01.00000.10001.09540.20001.18320.30001.26490.40001.34160.50001.41420.60001.48320.70001.54920.80001.61250.90001.67331.00001.7321步长h=0.2时，结果为:01.00000.20001.18320.40001.34170.60001.48330.80001.61251.00001.7321二、分析讨论由结果可知，改进的欧拉公式和四阶经典龙格库塔方法均能求解常微分方程的初值
6、问题，对于改进的欧拉格式，/?=0.1得到的结果比h=0.2更精确，同时,四阶经典龙格库塔方法也是如此，且对于实验所给的题目，/7=0.1时，四阶经典龙格库塔方法得到的结果和真实值最为接近。aa=1.7320508075688&plot(xl,y1,T,x2,y2,b)holdonplot(aa/mhr)legend('euler7nark47>F#确值Jgridon图一：h=0.2时，两种差分方法结果与真实值之间的比较步长为/?=0・2吋，两种常微分方程的差分方法结果比较同理可以得到步长为/心0.1时，两种常微分方程的差分方法结
7、果比较的图形。1.8图二：力=0.1时，两种差分方法结果与真实值之间的比较
