clear all
close all
clc
%比较Euler法、二阶泰勒级数法、Euler法的梯形算法与四阶RK算法-----
h=0.1; %step
x1=[0:h:1];
N=(max(x1)-min(x1))/h;
y0=1;  %the initial value
f=@(x,y)(y-2*x/y);
%--------------解析解------------------------
y_e=sqrt(1+2*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*(y_euler(n)-2*x1(n)/y_euler(n));
end
figure,subplot(2,1,1),plot(x1,y_e,'r',x1,y_euler,'b')
title('Euler法与精确结果比较图')
legend('精确解','Euler法结果')
subplot(2,1,2),plot(x1,abs(y_e-y_euler),'r')
title('Euler法与精确结果间的误差')
%----------二阶泰勒级数法----------
y_taylor=zeros(1,N+1);
y_taylor(1)=y0;
for n=1:N
    y_taylor(n+1)=y_taylor(n)+h*(y_taylor(n)-2*x1(n)/y_taylor(n))...
       +h^2/2*(-2/y_taylor(n)+(y_taylor(n)-2*x1(n)/y_taylor(n))*(1+2*x1(n)/y_taylor(n)^2));
end
%可见泰勒级数法比较复杂
figure,subplot(2,1,1),plot(x1,y_e,'r',x1,y_taylor,'b')
title('二阶泰勒级数法与精确结果比较图')
legend('精确解','二阶泰勒级数法结果')
subplot(2,1,2),plot(x1,abs(y_e-y_taylor),'r')
title('二阶泰勒级数法与精确结果间的误差')
%--------------改进的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,'b')
title('二阶RK算法与精确结果比较图')
legend('精确解','二阶RK算法结果')
subplot(2,1,2),plot(x1,abs(y_e-y_2rk),'r')
title('二阶RK算法与精确结果间的误差')
%--------------4th RK算法------------------------
y_4rk=zeros(1,N+1);
y_4rk(1)=y0;
for n=1:N
    k1=f(x1(n),y_4rk(n));
    k2=f(x1(n)+h/2,y_4rk(n)+h/2*k1);
    k3=f(x1(n)+h/2,y_4rk(n)+h/2*k2);
    k4=f(x1(n)+h,y_4rk(n)+h*k3);
    y_4rk(n+1)=y_4rk(n)+h/6*(k1+2*k2+2*k3+k4);
end  
figure,subplot(2,1,1),plot(x1,y_e,'r',x1,y_4rk,'b')
title('四阶RK算法与精确结果比较图')
legend('精确解','四阶RK算法结果')
subplot(2,1,2),plot(x1,abs(y_e-y_4rk),'r')
title('四阶RK算法与精确结果间的误差')