Euler法解微分方程-Matlab程序

合集下载

相关主题

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

%主程序main.m-----OK!

clear;

T=0.0;

Y=zeros(3,1);

Y(1)=1.0;Y(2)=1.0;Y(3)=1.0;

H1=0.05;M=3;EPS=1.0e-05;%EPS精度要求M方程个数H1拟定的输出步长

for i=1:10

[X,Y]=euler1(T,H1,Y,M,EPS)

T=T+H1;

End

%变步长euler方法

function [X,Y1]=euler1(T,H1,Y,M,EPS) %M-方程个数,EPS-精度,Y0-右端初值，T-自变量前一点值，H-步长

N=1;P=1+EPS;X=T;G=zeros(M,1);

H=H1;%H-在程序中要改变的步长H1－主程序中确定的输出步长

for i=1:M

C(i)=Y(i);

end

K1=zeros(M,1);K2=zeros(M,1);K3=zeros(M,1);K4=zeros(M,1);

while P>=EPS %变步长积分一步（H1）

for i=1:M

G(i)=Y(i);

Y(i)=C(i);

end

DT=H/N;

T=X;

%－－变步长积分过程

for j=1:N

K1=F(Y);

K2=F(Y+H/2*K1');

K3=F(Y+H/2*K2');

K4=F(Y+H*K3');

for i=1:M

Y(i)=Y(i)+H/6*(K1(i)+2*K2(i)+2*K3(i)+K4(i));

T=T+DT;

end

%－－－－－－－－－－－－－－－－－－－－－

P=0.0;

for i=1:M

Q=abs(Y(i)-G(i));

if Q>P

P=Q;

end

H=H/2.0;

N=N+N;

end

T=X;

X=T+H1;

Y1=Y;

%右端函数值function D=F(y)

D(1)=y(2);

D(2)=-1*y(1);

D(3)=y(3);