function fig7_8

%  comparison of various DE solvers
%  y' = f(t,y)  with   y(0) = y0 '

% logistic equation
t0=0;
tmax=3;
y0=0.01;
r=10;
a0=(1-y0)/y0;

% example 2
% t0=0;
% tmax=3;
% y0=1;
% exact=exp(tmax)*cos(tmax^2);

ii=0;
for im=1:5
    if im==1
        m=20;
    else
        m=10^im;
    end
    ii=ii+1; points(ii)=m;
    t=linspace(t0,tmax,m+1);
    k=t(2)-t(1);

    for i=1:m+1
        exact(i)=1/(1+a0*exp(-r*t(i)));
    end

    % euler
    y_euler=eulers(t,y0,k,m+1);
    e_euler(ii)=norm(exact-y_euler,inf);

    % backward euler
    y_beuler=beulers(t,y0,k,m+1);
    e_beuler(ii)=norm(exact-y_beuler,inf);

    % trap method
    y_trap=traps(t,y0,k,m+1);
    e_trap(ii)=norm(exact-y_trap,inf);

    %  RK2: Heun
    y_rk2H=rk2H(t,y0,k,m+1);
    e_rk2H(ii)=norm(exact-y_rk2H,inf);

    %  RK2: min error
    y_rk2E=rk2E(t,y0,k,m+1);
    e_rk2E(ii)=norm(exact-y_rk2E,inf);

end

clf
% get(gcf)
set(gcf,'Position', [25 1115 658 230])
subaxis(1,1,1,1,'MT',0.01,'MB',0.18,'MR',-0.02,'ML',0.05,'P',0.04)
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)

loglog(points,e_trap,'--bo','MarkerSize',8,'LineWidth',1.4)
hold on
set(gca,'ytick',[1e-9 1e-6 1e-3 1]);
axis([10 1e5 1e-9 1])
loglog(points,e_rk2H,'--kd','MarkerSize',8,'LineWidth',1.4)
loglog(points,e_rk2E,'--r*','MarkerSize',8,'LineWidth',1.4)
grid on
set(gca,'MinorGridLineStyle','none')

%  get(legend)
%  'Position',[left, bottom, width, height]
legend({' Trap',' Heun',' RK2-R'},'Location','SouthWest','FontSize',16,'FontWeight','bold')

xlabel('Number of Time Points')
ylabel('Error')

set(gca,'FontSize',16,'FontWeight','bold')
hold off

%exportgraphics(gcf,'/Users/mark/Desktop/rk2.eps')


% right-hand side of DE
function z=f(t,y)
r=10; z=r*y*(1-y);
%z=y-2*t*exp(t)*sin(t^2);

% euler method
function ypoints=eulers(t,y0,k,n)
y=y0;
ypoints=y0;
for i=2:n
    yy=y+k*f(t(i-1),y);
    ypoints=[ypoints, yy];
    y=yy;
end

% backward euler method
function ypoints=beulers(t,y0,h,n)
tol=5*eps;
y=y0; fold=f(t(1),y);
ypoints=y0;
for i=2:n
    %  secant method
    yb=y;         fb=yb-h*f(t(i),yb)-y;
    yc=y+2*h*(fold+f(t(i),y+h*fold));  fc=yc-h*f(t(i),yc)-y;
    while abs(yc-yb)>tol
        ya=yb; fa=fb;
        yb=yc; fb=fc;
        yc=yb-fb*(yb-ya)/(fb-fa);
        fc=yc-h*f(t(i),yc)-y;
    end
    y=yc; fold=f(t(i),y);
    ypoints=[ypoints, y];
end

% trap method
function ypoints=traps(t,y0,k,n)
tol=5*eps;
y=y0;
ypoints=y0;
fold=y+0.5*k*f(t(1),y);  % euler
yc=y0;  yb=1.001*y0;
for i=2:n
    %  secant method
    ya=1.01*yc;  fa=ya-0.5*k*f(t(i),ya)-fold;
    yb=yc;  fb=yb-0.5*k*f(t(i),yb)-fold;
    err=1;
    while err>tol
        yc=yb-fb*(yb-ya)/(fb-fa);
        fc=yc-0.5*k*f(t(i),yc)-fold;
        err=abs(yc-yb);
        ya=yb; fa=fb; yb=yc; fb=fc;
    end
    y=yc; fold=y+0.5*k*f(t(i),y);
    ypoints=[ypoints, y];
end

% RK2: Heun
function ypoints=rk2H(t,y0,h,n)
b=0.5;
a=1-b;  alpha=0.5/b;  beta=alpha;
y=y0;
ypoints=y0;
for i=2:n
    k1=h*f(t(i-1),y);
    k2=h*f(t(i-1)+alpha*h,y+beta*k1);
    yy=y+a*k1+b*k2;
    ypoints=[ypoints, yy];
    y=yy;
end

% RK2: min error
function ypoints=rk2E(t,y0,h,n)
b=3/4;
a=1-b;  alpha=2/3;  beta=alpha;
y=y0;
ypoints=y0;
for i=2:n
    k1=h*f(t(i-1),y);
    k2=h*f(t(i-1)+alpha*h,y+beta*k1);
    yy=y+a*k1+b*k2;
    ypoints=[ypoints, yy];
    y=yy;
end

% RK4 method
function ypoints=rk4s(t,y0,h,n)
y=y0;
ypoints=y0;
for i=2:n
    k1=h*f(t(i-1),y);
    k2=h*f(t(i-1)+0.5*h,y+0.5*k1);
    k3=h*f(t(i-1)+0.5*h,y+0.5*k2);
    k4=h*f(t(i),y+k3);
    yy=y+(k1+2*k2+2*k3+k4)/6;
    ypoints=[ypoints, yy];
    y=yy;
end




