clear all
close all
clc
%比较Euler法、Euler法的梯形算法
h=0.1;%step %改变步长，考察稳定性
x1=[0:h:0.5];
N=(max(x1)-min(x1))/h;
y0=1;  %the initial value
f=@(x,y)(-30*y);
%--------------解析解------------------------
y_e=exp(-30*x1);
%----------Euler算法------------------------------
y_euler=zeros(1,N+1);
y_euler(1)=y0;
for n=1:N %Euler法迭代公式
    y_euler(n+1)=y_euler(n)-h*30*y_euler(n);
end
figure,subplot(2,1,1),plot(x1,y_e,'r',x1,y_euler,'-ob')
title('Euler法与精确结果比较图')
legend('精确解','Euler法结果')
subplot(2,1,2),plot(x1,abs(y_e-y_euler),'b')
title('Euler法与精确结果误差曲线')
y_euler

%--------------改进的Euler法,二阶RK算法----------------
y_2rk=zeros(1,N+1);
y_2rk(1)=y0;
for n=1:N
    k1=f(x1(n),y_2rk(n));
    k2=f(x1(n)+h,y_2rk(n)+h*k1);
    y_2rk(n+1)=y_2rk(n)+h/2*(k1+k2);
end
figure,subplot(2,1,1),plot(x1,y_e,'r',x1,y_2rk,'-ob')
title('二阶RK算法与精确结果比较图')
legend('精确解','二阶RK算法结果')
subplot(2,1,2),plot(x1,abs(y_e-y_2rk),'b')
title('二阶RK算法与精确结果误差曲线')
y_2rk
