Euler and Runge Kutta methods

向前Euler方法

function Forward_Euler(fun,u_0,a,b,N)
%% 用 Forward-Euler法来求区间[a,b]的常微分数值解,N为区间[a,b]分割的份数
%u'(t)=f(t,u);
%u(a)=u_0;
%迭代格式 U_k+1=U_k+delta_t*f(t_k,U_k);
%显式格式


% step1 将区间[a,b]等分为N份
delta_t=(b-a)/N;
y=zeros(N+1,1);
b=zeros(N+1,1); % 初始化向量
%step2 迭代求解
for i=1:N+1
    b(i)=a+(i-1)*delta_t;
end
y(1)=u_0;
for i=1:N
    y(i+1)=y(i)+delta_t*feval(fun,b(i),y(i));
end
plot(b,y);


end
function f=fun(t,u)
f=u+t;
end

 

 

向后Euler 方法

function Backward_Euler(fun,u_0,a,b,N)
%% 用 Backward-Euler法来求区间[a,b]的常微分数值解,N为区间[a,b]分割的份数
%u'(t)=f(t,u);
%u(a)=u_0;
%迭代格式 U_k+1=U_k+delta_t*f(t_k+1,U_k+1);
% 隐式格式


% step1 将区间[a,b]等分为N份

delta_t=(b-a)/N;
y=zeros(N+1,1);
b=zeros(N+1,1); % 初始化向量

%step2 迭代求解
for i=1:N+1
    b(i)=a+(i-1)*delta_t;
end
y(1)=u_0;
for i=1:N
    %  iteration solve  y(i+1)=y(i)+delta_t*feval(fun,b(i+1),y(i+1));
    eps=10^(-2); %解的精度
    m1=0;temp=0;
    m2=1;
    while (abs(m1-m2)>eps)
        m1=temp;
        m2=y(i)+delta_t*feval(fun,b(i+1),m1);
        temp=m2;
    end
    y(i+1)=m2;
end
plot(b,y,'y');


end
function f=fun(t,u)
f=u+t;
end

梯形法

function trapezoid(fun,u_0,a,b,N)
%% 用 trapezoid法来求区间[a,b]的常微分数值解,N为区间[a,b]分割的份数
%u'(t)=f(t,u);
%u(a)=u_0;
%迭代格式 U_k+1=U_k+delta_t/2*(f(t_k,U_k)+f(t_k+1,u_k+1));
%显式格式


% step1 将区间[a,b]等分为N份
delta_t=(b-a)/N;
y=zeros(N+1,1);
b=zeros(N+1,1); % 初始化向量
%step2 迭代求解
for i=1:N+1
    b(i)=a+(i-1)*delta_t;
end
y(1)=u_0;
for i=1:N
    %  iteration solve  y(i+1)=y(i)+delta_t*feval(fun,b(i+1),y(i+1));
    eps=10^(-2); %解的精度
    m1=0;temp=0;
    m2=1;
    while (abs(m1-m2)>eps)
        m1=temp;
        m2=y(i)+delta_t/2*(feval(fun,b(i+1),m1)+feval(fun,b(i),y(i)));
        temp=m2;
    end
    y(i+1)=m2;
end
plot(b,y);


end
function f=fun(t,u)
f=u+t;
end



%---------------------example--------------------------
%trapezoid('fun',1,0,2,1000)

 

改进Euler 方法

function Improved_Euler(fun,u_0,a,b,N)
%% 用 Improved-Euler法来求区间[a,b]的常微分数值解,N为区间[a,b]分割的份数
%u'(t)=f(t,u);
%u(a)=u_0;
%迭代格式 U_k+1=U_k+delta_t*f(t_k,U_k);
%         U_k+1=U_k+delta_t/2*(f(t_k,U_k)+f(t_k+1,U_k+1));
% 显式格式  second order Runge-Kutta method


% step1 将区间[a,b]等分为N份
delta_t=(b-a)/N;
y=zeros(N+1,1);
b=zeros(N+1,1); % 初始化向量
%step2 迭代求解
for i=1:N+1
    b(i)=a+(i-1)*delta_t;
end
y(1)=u_0;
for i=1:N
    y(i+1)=y(i)+delta_t*feval(fun,b(i),y(i));
    y(i+1)=y(i)+delta_t*(feval(fun,b(i),y(i))+feval(fun,b(i+1),y(i+1)))/2;
end
plot(b,y,'r');

end
function f=fun(t,u)
f=t+u;
end

 

四阶Runge-Kutta方法

function RK4(fun,u_0,a,b,N)
%% 用 Fifth Order Runge Kutta来求区间[a,b]的常微分数值解,N为区间[a,b]分割的份数
%u'(t)=f(t,u);
%u(a)=u_0;
%迭代格式 U_k+1=U_k+delta_t/6*(K_1+2K_2+2K_3+K_4);
% K_1=f(t_k,U_k);
% K_2=f(t_k+delta_t*1/2，U_k+delta_t*K_1*1/2);
% K_3=f(t_k+delta_t*1/2，U_k+delta_t*K_2*1/2);
% K_4=f(t_k+delta_t,U_k+delta_t*K_3);
% 显式格式


% step1 将区间[a,b]等分为N份
delta_t=(b-a)/N;
y=zeros(N+1,1);
b=zeros(N+1,1); % 初始化向量
%step2 迭代求解
for i=1:N+1
    b(i)=a+(i-1)*delta_t;
end
y(1)=u_0;
for i=1:N
    K_1=feval(fun,b(i),y(i));
    K_2=feval(fun,b(i)+delta_t/2,y(i)+delta_t*K_1*1/2);
    K_3=feval(fun,b(i)+delta_t/2,y(i)+delta_t*K_2*1/2);
    K_4=feval(fun,b(i)+delta_t,y(i)+delta_t*K_3);
    y(i+1)=y(i)+delta_t/6*(K_1+2*K_2+2*K_3+K_4);
end
plot(b,y,'r');

end
function f=fun(t,u)
f=u*sin(u)+t*cos(u);
end