%% MCHI, cv6, file 1 (cv6_1) - Numerical solution to sode (A-B method)
%
% Adams-Bashforth method for solving a system of ordinary differential
% equations. It is a multistep solver. This exaple consits of a fourth
% order variation of this method. Thus first four steps have to be prepared
% by a different solver. To this I propose to use the already programmed
% fourth order Runge-Kutta method (see file cv4_1.m)
%
% This way we will be able to maintain the overall order of method (four)
%
% Please note, that being an explicit solver, A-B method is not A-stable
% and thus not suitable for solving stiff systems
%
% Author:       Martin Isoz
% Organisation: ICT Prague
% Date:         19. 10. 2014
%
% License: This code is published under MIT License, please do not abuse
% it.
%

% Note: dxdt = F(x) ... system of n ordinary differential equations

%% Caller 'script'
function caller
close all

% solving problem #1
Y0      = 1;
TSPAN   = [0 1];
H       = 0.0005;

% comparison of different methods
tstart = tic;
[TOUT,YOUTE] = eulerM(@problem1,TSPAN,Y0,H);                                %Euler - quick/innacurrate
timeE = toc(tstart);
tstart = tic;
[~,YOUTAB4]  =   ab4M(@problem1,TSPAN,Y0,H);                                %AB4(!!) - slow/acurrate
timeAB4 = toc(tstart);
tstart = tic;
[~,YOUTAB4O]  =   ab4MO(@problem1,TSPAN,Y0,H);                              %AB4 - quicker/acurrate
timeAB4O = toc(tstart);
tstart = tic;
[~,YOUTRK4]  =   rk4M(@problem1,TSPAN,Y0,H);                                %RK4 - slow/acurrate
timeRK4 = toc(tstart);

figure('Units','Pixels','Position',[100 50 1000 700],'Visible','on');       %open a figure window
subplot(211)                                                                %plot solution
    plot(TOUT,exp(TOUT),'k-');                                              %ploting the exact solution
    hold on;
    plot(TOUT,YOUTE,'ko');plot(TOUT,YOUTAB4,'ks');plot(TOUT,YOUTAB4O,'k^');
    plot(TOUT,YOUTRK4,'kd');                                                %ploting numerical approximations
    axis tight
    title('Solution of problem 1','FontSize',18,'FontWeight','Bold')
    xlabel('t');ylabel('y(t)');
    legend(...
        'Exact',...
        ['Euler t = ' num2str(timeE) ' s'],...
        [' AB4  t = ' num2str(timeAB4) ' s'],...
        [' AB4O t = ' num2str(timeAB4O) ' s'],...
        [' RK4  t = ' num2str(timeRK4) ' s'],'Location','Best');
subplot(212)                                                                %plot relative errors
    semilogy(TOUT,abs(YOUTE - exp(TOUT))./exp(TOUT).*100,'ko');hold on;     %euler method relative error
    semilogy(TOUT,abs(YOUTAB4 - exp(TOUT))./exp(TOUT).*100,'ks');           %  ab4 method relative error
    semilogy(TOUT,abs(YOUTRK4 - exp(TOUT))./exp(TOUT).*100,'kd');           %  rk4 method relative error
    axis tight
    title('Relative errors for problem 1','FontSize',18,'FontWeight','Bold')
    xlabel('t');ylabel('\Delta_R^%');
    legend(...
        'Euler, O(h^1)',...
        'Adams-Bashforth, O(h^4)',...
        'Runge-Kutta, O(h^4)','Location','Best');

% solving problem #2
Y0      = [2/3 1];
TSPAN   = [0 1];
H       = 0.1;

[TOUT,YOUTE] = eulerM(@problem2,TSPAN,Y0,H);
[~,YOUTAB4]  =  ab4MO(@problem2,TSPAN,Y0,H);

figure
plot(TOUT,exactSol2(TOUT));hold on;plot(TOUT,YOUTE,'o');plot(TOUT,YOUTAB4,'s');
legend('xExact','yExact','xEuler','yEuler','xAB4','yAB4','Location','Best')
axis tight
title('Solution of problem 2','FontSize',18,'FontWeight','Bold')

% solving problem #3
% this problem will be solved only using AB4M solver but for the different
% set of initial conditions

% create a matrix for different initial conditions
Y0M = [0.2:0.1:1;0.2:0.1:1]';                                                 %create matrix 2xN than transpose it to Nx2

% prepare the rest of inputs for solver
TSPAN   = [0 10];
H       = 0.1;

% prepare parameters for the model
alpha = 1;beta = 1;gamma = 1;delta = 1;
pars  = [alpha beta gamma delta];                                           %same structure as in problem3 function

% prepare figure window
figure;hold on

% solve the problem
for i = 1:numel(Y0M)/2
    Y0 = Y0M(i,:);                                                          %I will start in point (0.2,0.2),(0.4,0.4)...
    [~,YOUT]  =  ab4MO(@(t,x) problem3(t,x,pars),TSPAN,Y0,H);               %this time I am not interested in dependent variable
    plot(YOUT(:,1),YOUT(:,2),'.','Color',rand(1,3));                        %plot trajectory in the phase space
end

% format the graphics
title('Trajectories of problem3 in plane x_1-x_2','FontSize',18,'FontWeight','Bold')
xlabel('x_1');ylabel('x_2');

end

%% Methods

% Eulers method
function [TOUT,YOUT] = eulerM(ODEFUN,TSPAN,Y0,H)                            %with tspan we input time step (not adaptive)

% calculate needed # of steps
nSteps = norm(diff(TSPAN))./H+1;                                            %I will not end exactly in TEND

% allocate variables
TOUT = zeros(nSteps,1); YOUT = zeros(nSteps,numel(Y0));

% save initial condition
TOUT(1) = min(TSPAN);YOUT(1,:) = reshape(Y0,1,numel(Y0));

for i = 2:nSteps
    TOUT(i) = TOUT(i-1) + H;
    YOUT(i,:) = YOUT(i-1,:) + H*reshape(feval(ODEFUN,TOUT(i-1),YOUT(i-1,:)),1,numel(Y0));
end

end

% Runge-Kutta method (4th order) 
function [TOUT,YOUT] = rk4M(ODEFUN,TSPAN,Y0,H)                              %with tspan we input time step (not adaptive)
% calculate needed # of steps
nSteps = norm(diff(TSPAN))./H+1;                                            %I will not end exactly in TEND

% allocate variables
TOUT = zeros(nSteps,1); YOUT = zeros(nSteps,numel(Y0));

% save initial condition
TOUT(1) = min(TSPAN);YOUT(1,:) = reshape(Y0,1,numel(Y0));

for i = 2:nSteps
   TOUT(i) = TOUT(i-1) + H; 
   k1 = H*reshape(feval(ODEFUN,TOUT(i-1),YOUT(i-1,:)),1,numel(Y0));
   k2 = H*reshape(feval(ODEFUN,TOUT(i-1)+0.5*H,YOUT(i-1,:)+0.5*k1),1,numel(Y0));
   k3 = H*reshape(feval(ODEFUN,TOUT(i-1)+0.5*H,YOUT(i-1,:)+0.5*k2),1,numel(Y0));
   k4 = H*reshape(feval(ODEFUN,TOUT(i-1)+1*H,YOUT(i-1,:)+1*k3),1,numel(Y0));
   YOUT(i,:) = YOUT(i-1,:) + 1/6*(k1+k4) + 1/3*(k2+k3);
end
end

% Adams-Bashforts method (4th order)
function [TOUT,YOUT] = ab4M(ODEFUN,TSPAN,Y0,H)                              %with tspan we input time step (not adaptive)
% calculate needed # of steps
nSteps = norm(diff(TSPAN))./H+1;                                            %I will not end exactly in TEND

% allocate variables
TOUT = zeros(nSteps,1); YOUT = zeros(nSteps,numel(Y0));

% save initial condition
TOUT(1) = min(TSPAN);YOUT(1,:) = reshape(Y0,1,numel(Y0));

% method run starting - use 4th order RK method to get first 4 points of
% the solution
for i = 2:5
   TOUT(i) = TOUT(i-1) + H; 
   k1 = H*reshape(feval(ODEFUN,TOUT(i-1),YOUT(i-1,:)),1,numel(Y0));
   k2 = H*reshape(feval(ODEFUN,TOUT(i-1)+0.5*H,YOUT(i-1,:)+0.5*k1),1,numel(Y0));
   k3 = H*reshape(feval(ODEFUN,TOUT(i-1)+0.5*H,YOUT(i-1,:)+0.5*k2),1,numel(Y0));
   k4 = H*reshape(feval(ODEFUN,TOUT(i-1)+1*H,YOUT(i-1,:)+1*k3),1,numel(Y0));
   YOUT(i,:) = YOUT(i-1,:) + 1/6*(k1+k4) + 1/3*(k2+k3);
end

% the method itself - we use the previously obtained solution points to get
% the new ones
% y(n+1) = y(n) + h*sum_{i = 1}^p \beta_{qi}f_{n-i} (p-step method)
for i = 6:nSteps
    TOUT(i) = TOUT(i-1) + H;
    YOUT(i,:) = YOUT(i-1,:) + H/720*(...
        1901*reshape(feval(ODEFUN,TOUT(i-1),YOUT(i-1,:)),1,numel(Y0))...
        -2774*reshape(feval(ODEFUN,TOUT(i-2),YOUT(i-2,:)),1,numel(Y0)) + ...
        2616*reshape(feval(ODEFUN,TOUT(i-3),YOUT(i-3,:)),1,numel(Y0))...
        -1274*reshape(feval(ODEFUN,TOUT(i-4),YOUT(i-4,:)),1,numel(Y0)) +...
        251*reshape(feval(ODEFUN,TOUT(i-5),YOUT(i-5,:)),1,numel(Y0)));
end

end

% Adams-Bashforts method (4th order) - OPTIMIZED (at least somehow)
% I can store all the previously calculated RHS values in an auxiliary
% variable - AUX
function [TOUT,YOUT] = ab4MO(ODEFUN,TSPAN,Y0,H)                             %with tspan we input time step (not adaptive)
% calculate needed # of steps
nSteps = norm(diff(TSPAN))./H+1;                                            %I will not end exactly in TEND
nStart = 6;                                                                 %number of points needed to start the method

% allocate variables
TOUT = zeros(nSteps,1); YOUT = zeros(nSteps,numel(Y0));
AUX  = zeros(5,numel(Y0));                                                  %matrix to keep the known values in

% save initial condition
TOUT(1) = min(TSPAN);YOUT(1,:) = reshape(Y0,1,numel(Y0));

% method run starting - use 4th order RK method to get first 4 points of
% the solution
for i = 2:nStart
    TOUT(i) = TOUT(i-1) + H; 
   k1 = H*reshape(feval(ODEFUN,TOUT(i-1),YOUT(i-1,:)),1,numel(Y0));
   k2 = H*reshape(feval(ODEFUN,TOUT(i-1)+0.5*H,YOUT(i-1,:)+0.5*k1),1,numel(Y0));
   k3 = H*reshape(feval(ODEFUN,TOUT(i-1)+0.5*H,YOUT(i-1,:)+0.5*k2),1,numel(Y0));
   k4 = H*reshape(feval(ODEFUN,TOUT(i-1)+1*H,YOUT(i-1,:)+1*k3),1,numel(Y0));
   YOUT(i,:) = YOUT(i-1,:) + 1/6*(k1+k4) + 1/3*(k2+k3);
   AUX(i-1,:)  = k1;                                                        %fill in the auxiliary variable
end

% the method itself - we use the previously obtained solution points to get
% the new ones
% y(n+1) = y(n) + h*sum_{i = 1}^p \beta_{qi}f_{n-i} (p-step method)
for i = nStart+1:nSteps
    % update the auxiliary matrix
    AUX(1:4,:)= AUX(2:5,:);                                                 %create shape for the current values
    AUX(5,:)  = H*reshape(feval(ODEFUN,TOUT(i-1),YOUT(i-1,:)),1,numel(Y0)); %fill in current values of RHS
    % calculate the next step
    TOUT(i) = TOUT(i-1) + H;
    YOUT(i,:) = YOUT(i-1,:) + (...
        1901*AUX(5,:)...
        -2774*AUX(4,:) + ...
        2616*AUX(3,:)...
        -1274*AUX(2,:) +...
        251*AUX(1,:))./720;
end

end

%% Problem definitions
% first equation to be solved (x in R)
function dxdt = problem1(~,x)                                               %I do not need to input t -> ~
dxdt = x;
end

% second equation to be solved (x in R2)
function dxdt = problem2(~,x)                                               %%I do not need to input t -> ~

dxdt = zeros(numel(x),1);                                                   %allocate space for output (general)

dxdt(1) = x(1) + x(2).^2;                                                   %problem specific
dxdt(2) = -x(2);                                                            %problem specific
end

% exact solution to the problem 2
function xExact = exactSol2(t)

xExact = zeros(numel(t),2);
xExact(:,1) = exp(t) - exp(-2*t)./3;
xExact(:,2) = exp(-t);
end

function dxdt = problem3(~,x,pars)
% Note: This is a SODE defining the classical predator-prey model
% (Lotka-Voltera equations)

dxdt = zeros(numel(x),1);

% unwrap parameters
alpha = pars(1);
beta  = pars(2);
gamma = pars(3);
delta = pars(4);

% calculate the derivations
dxdt(1) = x(1)*(alpha - beta*x(2));
dxdt(2) = -x(2).*(gamma - delta*x(1));
end