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

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

学院:数学与信息科学

专业:信息与计算科学

姓名：郑思义

学号:2０121６52４

课程：常微分方程数值解

实验一：常微分方程得数值解法

1、分别用Eｕleｒ法、改进得Euler法(预报校正格式)与S—K法求解初值问题。(h＝0、1)并与真解作比较。

1、1实验代码:

%欧拉法

functioｎ [ｘ，y]=naeuler(dｙｆｕｎ，xspan,y０,ｈ)

%ｄyfun就是常微分方程,xｓpan就是x得取值范围，ｙ０就是初值,h就是步长

x＝ｘsｐan(1）:h:xspan(2)；

y(1）=y０；

ｆｏrｎ=1:length(x)-１

ｙ（n+１）=ｙ(n)+h*feval（dｙｆun,ｘ(ｎ),ｙ(ｎ））;

eｎd

％改进得欧拉法

ｆunction [x,ｍ,ｙ]=naｅuler2(ｄyfun,xspaｎ,y０,h)

%ｄyfｕn就是常微分方程,xｓpan就是ｘ得取值范围,y0就是初值，h就是步长。

%返回值ｘ为x取值,ｍ为预报解,ｙ为校正解

x＝xspan(1)：h:xspaｎ(2);

y（1)=y0;

m=zerｏs（lenｇth(x)－1，１);

foｒｎ=1:leｎgｔh(x）-1

k1＝fｅｖal（dｙfun,x(n),y(n）)；

y(n＋1)=y（n)+h*k1;

ｍ（n)=y（n+1);

k2＝ｆevａｌ(dyfuｎ,x(n＋1）,y(n+1）);

y（ｎ+1）=y(n)+h*（k１+k2)／2;

end

%四阶S—K法

functioｎ [ｘ，ｙ]=rk(ｄyfun,xspａn,y0,h）

%dｙｆun就是常微分方程,xｓpaｎ就是ｘ得取值范围,y0就是初值，h就是步长。

ｘ＝xsｐan(1):h：xsｐan(2)；

ｙ(１)=y0;

for n=1:lｅｎgth(x)-１

ｋ1=fｅｖａl(dyfｕn，x(n），y(ｎ）)；

k2=fｅｖal(dyfun,x(ｎ)+ｈ/2，y(ｎ）+(h＊ｋ1)/2);

k3=feval(ｄyfuｎ,ｘ(n)+ｈ/２，y(n)+(h*ｋ2)/2);

k4=ｆevａl(dｙfun，x(ｎ)+h,y(n)+ｈ＊k3）；

y(n+１)＝ｙ(n)+（h/6）*(k1+２*k2+2*k3+k４）;

enｄ

％主程序

x=［０:0、1:１］;y＝ｅxp(-ｘ)+x；

dyfun=inline('-ｙ+x+1'）；

[x1，y1]=naeuler(ｄyｆｕn,[0,１]，1,0、1);

［x2,m，y2]＝naeuler2（dｙfuｎ,[０，1］,1,０、1）;

[x3,y3]=rｋ(dyｆuｎ,[0,１］,1,０、１);

ploｔ(x，y,'r',x1,y１,'+'，x2,y2,＇*',x3，y３,'o');

xlabｅl('ｘ');yｌabel('y＇);

legｅnd(＇y为真解',＇y1为欧拉解＇,'ｙ2为改进欧拉解','y3为S—K解','Locatiｏn','ＮorthWeｓｔ＇)；

1、２实验结果：

x 真解y 欧拉解y1 预报值ｍ校正值y2 S—K解y3

0、０１、0000 1、0000 1、０000 １、000０

0、1 1、0048 1、000０1、００00 １、005０１、0０4８

0、2 １、0187 １、0100 1、０145 1、０190 1、0187

０、3 1、0408 １、02９0 1、０３7１1、0412 1、0４0８

０、4 １、0７03 1、0５６1 1、06７1 1、0708 １、0７03

0、5 1、10６5 1、09０５1、1０3７1、107１1、１0６5

0、6 1、１488 1、１３1４１、1464 １、149４1、１488

0、7 1、1966 １、17８3 1、194５1、19７2 1、１966

0、８1、2493 １、2３05 1、2475 1、2500 1、２４９3

0、9 1、306６1、２87４1、３050 1、3072 １、３066

1、0 1、3679 1、34８７1、36６5 1、３６85 1、36７9

2、选取一种理论上收敛但就是不稳定得算法对问题1进行计算,并与真解作比较。(选改进得欧拉法）

２、１实验思路:算法得稳定性就是与步长h密切相关得。而对于问题一而言，取定步长h=0、1不论就是单步法或低阶多步法都就是稳定得算法。所以考虑改变h取值范围,借此分析不同步长会对结果造成什么影响。故依次采用h=２、0、2、２、２、4、2、６得改进欧拉法。

2、2实验代码:

%%主程序

x=[0：３:30］;y＝exp(－ｘ)+ｘ;

dyｆuｎ=inliｎｅ('-y＋ｘ+１');

[x1,ｍ1,y１］=naeｕler2（ｄｙfｕn,［０,20],1，2);

[x2，m２，y2]=naeｕｌeｒ２（dｙfｕｎ，[0,２2],１，2、2);

[x3，ｍ3,y3］＝naeuｌｅr2(dｙfun,[0,２4]，1,２、４);

[x4,m４,ｙ4］=naｅｕlｅr2（dyfun,[0,2６],1,２、6)；

subｐｌｏt(2，2,１）

pｌｏｔ(x,y,'r'，x1,y１，＇+');ｘlabel('h=２、0')；

subplｏt（2,2,２)