资源描述：
《常微分方程数值解实验报告材料》由会员上传分享，免费在线阅读，更多相关内容在应用文档-天天文库。
1、实用标准文案常微分方程数值解实验报告学院：数学与信息科学专业：信息与计算科学姓名：郑思义学号：201216524课程：常微分方程数值解文档实用标准文案实验一：常微分方程的数值解法1、分别用Euler法、改进的Euler法（预报校正格式）和S—K法求解初值问题。（h=0.1）并与真解作比较。1.1实验代码：%欧拉法function[x,y]=naeuler(dyfun,xspan,y0,h)%dyfun是常微分方程，xspan是x的取值范围，y0是初值，h是步长x=xspan(1):h:xspan(2);y(1)=y0;forn=1:
2、length(x)-1y(n+1)=y(n)+h*feval(dyfun,x(n),y(n));end%改进的欧拉法function[x,m,y]=naeuler2(dyfun,xspan,y0,h)%dyfun是常微分方程，xspan是x的取值范围，y0是初值，h是步长。%返回值x为x取值，m为预报解，y为校正解x=xspan(1):h:xspan(2);y(1)=y0;m=zeros(length(x)-1,1);forn=1:length(x)-1k1=feval(dyfun,x(n),y(n));y(n+1)=y(n)+h*
3、k1;m(n)=y(n+1);k2=feval(dyfun,x(n+1),y(n+1));y(n+1)=y(n)+h*(k1+k2)/2;end%四阶S—K法function[x,y]=rk(dyfun,xspan,y0,h)%dyfun是常微分方程，xspan是x的取值范围，y0是初值，h是步长。x=xspan(1):h:xspan(2);y(1)=y0;forn=1:length(x)-1k1=feval(dyfun,x(n),y(n));k2=feval(dyfun,x(n)+h/2,y(n)+(h*k1)/2);k3=fev
4、al(dyfun,x(n)+h/2,y(n)+(h*k2)/2);k4=feval(dyfun,x(n)+h,y(n)+h*k3);文档实用标准文案y(n+1)=y(n)+(h/6)*(k1+2*k2+2*k3+k4);end%主程序x=[0:0.1:1];y=exp(-x)+x;dyfun=inline('-y+x+1');[x1,y1]=naeuler(dyfun,[0,1],1,0.1);[x2,m,y2]=naeuler2(dyfun,[0,1],1,0.1);[x3,y3]=rk(dyfun,[0,1],1,0.1);pl
5、ot(x,y,'r',x1,y1,'+',x2,y2,'*',x3,y3,'o');xlabel('x');ylabel('y');legend('y为真解','y1为欧拉解','y2为改进欧拉解','y3为S—K解','Location','NorthWest');1.2实验结果：x真解y欧拉解y1预报值m校正值y2S—K解y30.01.00001.00001.00001.00000.11.00481.00001.00001.00501.00480.21.01871.01001.01451.01901.01870.31.04081
6、.02901.03711.04121.04080.41.07031.05611.06711.07081.07030.51.10651.09051.10371.10711.10650.61.14881.13141.14641.14941.14880.71.19661.17831.19451.19721.19660.81.24931.23051.24751.25001.24930.91.30661.28741.30501.30721.30661.01.36791.34871.36651.36851.3679文档实用标准文案1、选取一种理
7、论上收敛但是不稳定的算法对问题1进行计算，并与真解作比较。（选改进的欧拉法）2.1实验思路：算法的稳定性是与步长h密切相关的。而对于问题一而言，取定步长h=0.1不论是单步法或低阶多步法都是稳定的算法。所以考虑改变h取值范围，借此分析不同步长会对结果造成什么影响。故依次采用h=2.0、2.2、2.4、2.6的改进欧拉法。2.2实验代码：%%主程序x=[0:3:30];y=exp(-x)+x;dyfun=inline('-y+x+1');[x1,m1,y1]=naeuler2(dyfun,[0,20],1,2);[x2,m2,y2]=
8、naeuler2(dyfun,[0,22],1,2.2);[x3,m3,y3]=naeuler2(dyfun,[0,24],1,2.4);[x4,m4,y4]=naeuler2(dyfun,[0,26],1,2.6);subplot(
