微分方程求解

1. 欧拉方法

经过计算：

　　

求[a,b]上的初值问题 y'=f(t,y)，y(a)=y₀的近似解。

function E=euler(f,a,b,ya,M)
%Input  - f is the function entered  as a string 'f'
%       - a and b are the left and right end points
%       - ya is the initial condition y(a)
%       - M is the number of steps
%Output - E=[T' Y'] where T is the vector of abscissas and
%         Y is the vector of ordinates
h=(b-a)/M;
Y=zeros(1,M+1);
T=a:h:b;
Y(1)=ya;
for j=1:M
    Y(j+1)=Y(j+h*feval(f,T(j),Y(j)));
end
E=[T' Y'];
end

 2.四阶龙格-库塔法

计算公式为：

　　

求[a,b]上的初值问题 y'=f(t,y)，y(a)=y₀的近似解。

function R=rk4(f,a,b,ya,M)
%Input  - f is the function entered  as a string 'f'
%       - a and b are the left and right end points
%       - ya is the initial condition y(a)
%       - M is the number of steps
%Output - R=[T' Y'] where T is the vector of abscissas and
%         Y is the vector of ordinates
%
h=(b-a)/M;
T=zeros(1,M+1);
U=zeros(4,M+1);
T=a:h:b; 
Y(1)=ya;
for j=1:M
    k1=h*feval(f,T(j),Y(j));
    k2=h*feval(f,T(j)+h/2,Y(j)+k1/2);
    k3=h*feval(f,T(j)+h/2,Y(j)+k2/2);
    k4=h*feval(f,T(j)+h,Y(j)+k3);
    Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)/6;
end
R=[T' Y'];
end