%% MCHI, cv6, file 2  (cv6_2) - Numerical solution to sode (P-C method)
%
% Explicit solvers for SODE are not A-stable and thus not well for stiff
% systems. One of the way to overcome this is to use some kind of
% correction for the value of x_n proposed by the explicit solver. This
% correction may be provided by embeded implicit solver.
%
% These methods are called predictor-corrector. For details on them see the
% lecture notes from numerical methods (Kubicek, Dubcova..., Numericke
% metody a algoritmy, VSCHT Praha)
% 
% In this script, wee will used a combination of Adams-Bashforth (explicit)
% and Adams-Moulton (implicit) solvers. We also will need to use the Newton
% methods for solving a system of NAE created by the implicit solver. For
% this we wil take the next step proposed by the explicit solver as an
% initial guess.
%
% To start the algorithm, we may either use RK4 method with a very short
% step and hope for the best or we can use any 1 step implicit solver.
%
%
% 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 5];
H       = 0.2;

% 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;
[~,YOUTPC4]  =  PC4M(@problem1,TSPAN,Y0,H);                              %AB4 - quicker/acurrate
timePC4 = 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,YOUTPC4,'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'],...
        [' PC4  t = ' num2str(timePC4) ' 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(YOUTPC4 - exp(TOUT))./exp(TOUT).*100,'k^');           %  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)',...
        'Pred.-Corrector, O(h^{??})',...
        'Runge-Kutta, O(h^4)','Location','Best');



% 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]  =   PC4M(@(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

% Reference 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

% Possible simplifications of known methods using submethods
% Runge-Kutta method (4th order) - with submethod predStart
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; 
   YOUT(i,:) = predStart(ODEFUN,TOUT(i-1),YOUT(i-1,:),H);
end
end

% Adams-Bashforth method (4th order) - with submethods predStart a predMain
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
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;
    [YOUT(i,:),AUX(i-1,:)] = predStart(ODEFUN,TOUT(i-1),YOUT(i-1,:),H);
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
    TOUT(i) = TOUT(i-1) + H;
    [YOUT(i,:),AUX] = predMain(ODEFUN,TOUT(i-1),YOUT(i-1,:),H,AUX);
end
end

% Predictor-Corrector method itself
% For this method we use:
% 1. submethod predStart, which is based on RK4 method. This method is used
% as a predictor method to start the main solver.
% 2. submethod predMain, which is based on AB4 method. This method is used
% as the predictor method for the main solver.
% 3. submethod corrMain, which is a simple iteration scheme. This method is 
% used as a corrector in each solver step
function [TOUT,YOUT] = PC4M(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;
    [YOUT(i,:),AUX(i-1,:)] = predStart(ODEFUN,TOUT(i-1),YOUT(i-1,:),H);          %get the guess for next step from predictor
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
    TOUT(i) = TOUT(i-1) + H;
    [YVOL,AUX] = predMain(ODEFUN,TOUT(i-1),YOUT(i-1,:),H,AUX);              %get the guess for next step from predictor
    YOUT(i,:)  = corrMain(ODEFUN,TOUT(i-1),YOUT(i-1,:),YVOL,H,AUX,5);       %refine via corrector
end
end

% Used sub-methods

% Predictor for starting the algorithm - predStart
%
% Runge-Kutta method (4th order) - for starting the algoritm - this method
% is used as a predictor for starting the calculation - I need to output
% also the value of k1 to fill the auxiliary matrix for the Adams-Bashforth
% method used as the predictor during a main body of the calculation
function [YOUT, K1] = predStart(ODEFUN,TN,YN,H)                             %with tspan we input time step (not adaptive)
% This function returns only next value of Y, it is not a solver by itself,
% merely a calculation a one single next step based on the current solution
% (RK4 is a single step method)

% preparation - I only reshape the inptut the input variable
% I actually do not need this step as I am the only one controling the
% shape of YN (it is not accessible to any user-defined function)
YN = reshape(YN,1,numel(YN));

% calculation of the next step of the solution
K1 = H*reshape(feval(ODEFUN,TN,YN),1,numel(YN));                            %!!This is also an output variable
k2 = H*reshape(feval(ODEFUN,TN+0.5*H,YN+0.5*K1),1,numel(YN));
k3 = H*reshape(feval(ODEFUN,TN+0.5*H,YN+0.5*k2),1,numel(YN));
k4 = H*reshape(feval(ODEFUN,TN+1*H,YN+1*k3),1,numel(YN));
YOUT = YN + 1/6*(K1+k4) + 1/3*(k2+k3);
end

% Predictor for the main calculation itself - predMain
%
% Adams-Bashforth method of 4th order. 
% Once I have all the needed previous values 
% (see the matrix AUX in the file cv6_1.m), I can use them to speed up
% the calculation - I do not need to calculate yN+1 by RK4 method (which
% is quite expensive), but I can use the knowledge of previous points and
% use AM4 method which has (roughly) the same precision and is
% substantially cheaper

function [YOUT,AUX] = predMain(ODEFUN,TN,YN,H,AUX)
% This function returns only next value of Y, it is not a solver by itself,
% merely a calculation a one single next step based on the current solution
% and on the stored values of RHS vector field

% update the auxiliary matrix
    AUX(1:4,:)= AUX(2:5,:);                                                 %create shape for the current values
    AUX(5,:)  = H*reshape(feval(ODEFUN,TN,YN),1,numel(YN));                 %fill in current values of RHS
    % calculate the next step
    YOUT = YN + (1901*AUX(5,:) - 2774*AUX(4,:) + 2616*AUX(3,:)...
        - 1274*AUX(2,:) + 251*AUX(1,:))./720;
end

% Corrector for the main calculation itself - corrMain
function YOUT = corrMain(ODEFUN,TN,YN,YVOL,H,AUX,nIter)
% for specified number of iterations, nIter, this functions uses a simple
% iteration scheme (see page 69 in lecture notes Numericke metody a
% algoritmy) to refine the initial guess from the predictor method
%
% be aware, that this method converges only for small method steps (H)
for j = 1:nIter
    YVOL = YN + ...                                                         %previous step
        H*251/720*reshape(feval(ODEFUN,TN,YVOL),1,numel(YN)) +...           %current step
        (646*AUX(5,:) - 264*AUX(4,:) + ...                                  %previous steps
        106*AUX(3,:) - 19*AUX(2,:))./720;
end
YOUT = YVOL;                                                                %update current step
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

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