常微分方程数值解实验探索.doc

常微分方程数值解实验报告 学院：数学与信息科学 专业：信息与计算科学 姓名：郑思义 学号：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 校正值y2 S—K解y3 0.0 1.0000 1.0000 1.0000 1.0000 0.1 1.0048 1.0000 1.0000 1.0050 1.0048 0.2 1.0187 1.0100 1.0145 1.0190 1.0187 0.3 1.0408 1.0290 1.0371 1.0412 1.0408 0.4 1.0703 1.0561 1.0671 1.0708 1.0703 0.5 1.1065 1.0905 1.1037 1.1071 1.1065 0.6 1.1488 1.1314 1.1464 1.1494 1.1488 0.7 1.1966 1.1783 1.1945 1.1972 1.1966 0.8 1.2493 1.2305 1.2475 1.2500 1.2493 0.9 1.3066 1.2874 1.3050 1.3072 1.3066 1.0 1.3679 1.3487 1.3665 1.3685 1.3679 选取一种理论上收敛但是不稳定的算法对问题1进行计算，并与真解作比较。（选改进的欧拉法） 2.1实验思路：算法的稳定性是与步长h密切相关的。而对于问题一而言，取定步长h=0.1不论是单步法或低阶多步法都是稳定的算法。所以考虑改变h取值范围，借此分析不同步长会对结果造成什么影响。故依次