function fig7_9

%  comparison of the RK4 and other RK4s

global icase

clf
% get(gcf)
set(gcf,'Position', [9 870 742 475])
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)
subaxis(2,1,1,1,'MT',0.014,'MB',0.075,'MR',-0.01,'ML',0.06,'P',0.03)

% logistic eq
icase=1;
tmax=1;
y0=0.1;
a0=(1-y0)/y0;
r=10;
% exact=1/(1+a0*exp(-r*t));
M=[10 100 1000 1e4];
for im=1:length(M)
    m=M(im);
    points(im)=m;
    t=linspace(0,tmax,m+1);
    k=t(2)-t(1);

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

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

    % Labatto RK4
    y_rk4L=rk4L(t,y0,k,m+1);
    e_rk4L(im)=norm(exact-y_rk4L,inf);

    %  RK4
    y_rk4=rk4(t,y0,k,m+1);
    e_rk4(im)=norm(exact-y_rk4,inf);

    %  RK4 Ralston min error
    y_rk4R=rk4R(t,y0,k,m+1);
    e_rk4R(im)=norm(exact-y_rk4R,inf);

end

loglog(points,e_rk4,'--bo','MarkerSize',8,'LineWidth',1.4)
hold on
loglog(points,e_rk4L,'--r*','MarkerSize',8,'LineWidth',1.4)
loglog(points,e_trap,'--kd','MarkerSize',8,'LineWidth',1.4)
grid on
%axis([10 1e3 1e-8 1])
set(gca,'ytick',[1e-15 1e-12 1e-9 1e-6 1e-3 1])
legend({' RK4-C',' RK4-L',' Trap'},'Location','SouthWest','FontSize',16,'FontWeight','bold')

% loglog(points,e_rk4R,'--gs','MarkerSize',8,'LineWidth',1.4)

%xlabel('Number of Time Points')
ylabel('Error')
set(gca,'FontSize',16,'FontWeight','bold')
hold off


subaxis(2,1,1,2)
% forced problem
icase=2;
y0=1; tmax=3; z=5;
% exact=exp(t)*cos(z*t^2);

for im=1:length(M)
    m=M(im);
    points(im)=m;
    t=linspace(0,tmax,m+1);
    k=t(2)-t(1);

    for i=1:m+1
        exact2(i)=exp(t(i))*cos(z*t(i)^2);
    end

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

    % Labatto RK4
    y_rk4L2=rk4L(t,y0,k,m+1);
    e_rk4L2(im)=norm(exact2-y_rk4L2,inf);

    %  RK4
    y_rk42=rk4(t,y0,k,m+1);
    e_rk42(im)=norm(exact2-y_rk42,inf);

    %  RK4 Ralston min error
    y_rk4R2=rk4R(t,y0,k,m+1);
    e_rk4R2(im)=norm(exact2-y_rk4R2,inf);

end

loglog(points,e_rk42,'--bo','MarkerSize',8,'LineWidth',1.4)
hold on
loglog(points,e_rk4L2,'--r*','MarkerSize',8,'LineWidth',1.4)
loglog(points,e_trap,'--kd','MarkerSize',8,'LineWidth',1.4)
grid on
%axis([10 1e3 1e-8 1])
set(gca,'ytick',[1e-15 1e-12 1e-9 1e-6 1e-3 1])
legend({' RK4-C',' RK4-L',' Trap'},'Location','SouthWest','FontSize',16,'FontWeight','bold')

% loglog(points,e_rk4R2,'--gs','MarkerSize',8,'LineWidth',1.4)

xlabel('Number of Time Points')
ylabel('Error')
set(gca,'FontSize',16,'FontWeight','bold')
hold off

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


function g=f(t,y)
global icase
if icase==1
    g=10*y*(1-y);
else
    z=5; g=y-2*z*t*exp(t)*sin(z*t^2);
end


%  Labatto RK4
function ypoints=rk4L(t,y0,k,n)
a2=(1-1/sqrt(5))/2;
a3=(1+1/sqrt(5))/2;
b21=a2;
c1= 1/12;
c2= 5/12;
c3= 5/12;
c4= 1/12;
b31 = -1/4-(3/20)*sqrt(5);
b32 = (1/4)*sqrt(5)+3/4;
b41 = -1/4+(5/4)*sqrt(5);
b42 = -5/4-(3/4)*sqrt(5);
b43 = 5/2-(1/2)*sqrt(5);
ypoints=y0;
y=y0;
for i=2:n
    k1=k*f(t(i-1),y);
    k2=k*f(t(i-1)+a2*k,y+b21*k1);
    k3=k*f(t(i-1)+a3*k,y+b31*k1+b32*k2);
    k4=k*f(t(i),y+b41*k1+b42*k2+b43*k3);
    yy=y+c1*k1+c2*k2+c3*k3+c4*k4;
    ypoints=[ypoints, yy];
    y=yy;
end


% RK4 method
function ypoints=rk4(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

% 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);
yc=y0;  yb=1.0001*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

% RK4 min error (ralston)
function ypoints=rk4R(t,y0,h,n)
y=y0;
ypoints=y0;
a2=2/5;
a3=7/8-3*sqrt(5)/16;
b21=a2;
c1= 1/2 +  (1-2*(a2+a3))/(12*a2*a3);
c2= (2*a3-1)/(12*a2*(a3-a2)*(1-a2));
c3= (1-2*a2)/(12*a3*(a3-a2)*(1-a3));
c4= 1/2 +  (2*( a2+a3)-3)/(12*(1-a2)*(1-a3));
b32= a3*(a3-a2)/(2*a2*(1-2*a2));
b42= (1- a2)*( a2+ a3-1-(2*a3-1)^2)/(2*a2*( a3- a2)*( 6*a2*a3-4*( a2+ a3)+3));
b43= (1-2*a2)*(1- a2)*(1- a3)/( a3*( a3- a2)*( 6*a2*a3-4*( a2+ a3)+3));
b31 =a3- b32;
b41=1-b42- b43;
for i=2:n
    k1=h*f(t(i-1),y);
    k2=h*f(t(i-1)+a2*h,y+b21*k1);
    k3=h*f(t(i-1)+a3*h,y+b31*k1+b32*k2);
    k4=h*f(t(i),y+b41*k1+b42*k2+b43*k3);
    yy=y+c1*k1+c2*k2+c3*k3+c4*k4;
    ypoints=[ypoints, y];
    y=yy;
end
%e_rkR(im)=abs(exact-yy);




