欧拉法matlab程序学习课件.doc

合集下载

相关主题

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

1.Euler法

function [x,y]=naeuler(dyfun,xspan,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

x=x';y=y';

x1=0:0.2:1;y1=(1+2*x1).^0.5;

plot(x,y,x1,y1)

>> dyfun=inline('y-2*x/y');

[x,y]=naeuler(dyfun,[0,1],1,0.2);[x,y]

ans =

0 1.0000

0.2000 1.2000

0.4000 1.3733

0.6000 1.5315

0.8000 1.6811

1.0000 1.8269

2.隐式Euler法

function [x,y]=naeulerb(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2);

y(1)=y0;

for n=1:length(x)-1

y(n+1)=iter(dyfun,x(n+1),y(n),h);

end

x=x';y=y';

x1=0:0.2:1;y1=(1+2*x1).^0.5;

plot(x,y,x1,y1)

function y=iter(dyfun,x,y,h)

y0=y;e=1e-4;K=1e+4;

y=y+h*feval(dyfun,x,y);

y1=y+2*e;k=1;

while abs(y-y1)>e

y1=y;

y=y0+h*feval(dyfun,x,y);

k=k+1;

if k>K

error('迭代发散');

end

>> dyfun=inline('y-2*x/y');

[x,y]=naeulerb(dyfun,[0,1],1,0.2);[x,y]

ans =

0 1.0000

0.2000 1.1641

0.4000 1.3014

0.6000 1.4146

0.8000 1.5019

1.0000 1.5561

3.改进Euler法

function [x,y]=naeuler2(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2);

y(1)=y0;

for n=1:length(x)-1

k1=feval(dyfun,x(n),y(n));

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

k2=feval(dyfun,x(n+1),y(n+1));

y(n+1)=y(n)+h*(k1+k2)/2;

end

x=x';y=y';

x1=0:0.2:1;y1=(1+2*x1).^0.5;

plot(x,y,x1,y1)

>> dyfun=inline('y-2*x/y');

[x,y]=naeuler2(dyfun,[0,1],1,0.2);[x,y] ans =

0 1.0000

0.2000 1.1867

0.4000 1.3483

0.6000 1.4937

0.8000 1.6279

1.0000 1.7542