- 最简单求解一个微分方程数值解:Euler法
function [x,y]=Euler_method(dufun,span,h,x0,y0) %EuLer格式, %求解方程y‘=dufun(x,y);其中x \in[a,b];y0为初始值;n为自变量的离散个数;y为求解结果 x=span(1):h:span(2); n=length(x); y=zeros(1,n);%存放数值的解 x(1)=x0; y(1)=y0; for i=1:n-1 y(i+1)=y(i)+h.*feval(dufun,x(i),y(i));%Euler格式 end
2.改进的Euler格式
function[x,y]=Gaijin_Euler(func,span,h,x0,y0)
%改进的Euler格式;
%span为区间
%y(x0)=y0;
%n表示区间等份;
x=span(1):h:span(2);
n=length(x);
y=zeros(1,n);
x(1)=x0;
y(1)=y0;
for i=1:n
y1=y(i)+h*feval(func,x(i),y(i));
y2=y(i)+h*feval(func,x(i+1),y1);
y(i)=(y1+y2)/2;
end
方法二:常见的Rungerkutta3与Rungerkutta4格式
1.Rungerkutta3function [x,y]=Rungekutta3(f1,span,h,x0,y0) x=span(1):h:span(2); n=length(x); y=zeros(1,n); x(1)=x0; y(1)=y0; for j=1:n-1 K1=feval(f1,x(j),y(j)); K2=feval(f1,x(j)+h/2,y(j)+h/2*K1); K3=feval(f1,x(j)+h,y(j)-h*K1+h*2*K2); y(j+1)=y(j)+(h/6)*(K1+4*K2+K3); end
1.Rungerkutta4function [x,y]=Rungekutta4(f1,span,h,x0,y0) x=span(1):h:span(2); n=length(x); y=zeros(1,n); x(1)=x0; y(1)=y0; for j=1:n-1 K1=feval(f1,x(j),y(j)); K2=feval(f1,x(j)+h/2,y(j)+h/2*K1); K3=feval(f1,x(j)+h/2,y(j)+h/2*K2); K4=feval(f1,x(j)+h,y(j)+h*K3); y(j+1)=y(j)+(h/6)*(K1+2*K2+2*K3+K4); end
时间: 2024-10-16 22:00:08