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

常微分方程数值解实验报告 学院：数学与信息科学 专业：信息与计算科学 姓名：郑思义 学号：201216524 课程：常微分方程数值解 实验一：常微分方程的数值解法 分别用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; for n=1:length(x)-1 y(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); for n=1:length(x)-1 k1=feval(dyfun,x(n),y(n)); y(n+1)=y(n)+h*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; for n=1:length(x)-1 k1=feval(dyfun,x(n),y(n)); k2=feval(dyfun,x(n)+h/2,y(n)+(h*k1)/2); k3=feval(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); plot(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.0 1.00001.00001.00001.00000.1 1.00481.00001.00001.00501.00480.2 1.01871.01001.01451.01901.01870.3 1.04081.02901.03711.04121.04080.4 1.07031.05611.06711.07081.07030.5 1.10651.09051.10371.10711.10650.6 1.14881.13141.14641.14941.14880.7 1.19661.17831.19451.19721.19660.8 1.24931.23051.24751.25001.24930.9 1.30661.28741.30501.30721.30661.0 1.36791.34871.36651.36851.3679 选取一种理论上收敛但是不稳定的算法对问题1进行计算，并与真解作比较。（选改进的欧拉法） 2.1实验思路：算法的稳定性是与步长h密切相关的。而对于问题一而言，取定步长h=0.1不论是单步法或低阶多步法都是稳定的算法。所以考虑改变h取值范围，借此分析不同步长会对结果造成什么影响。故依次