% ode.m - Matthieu Aussal - Copyright 2018 - Licence CC BY-NC-SA 4.0
% CMAP - Ecole Polytechnique, Route de Saclay, 91128 Palaiseau CEDEX France
% For any informations, please contact  : matthieu.aussal@polytechnique.edu

clear all
close all
clc

%%%%%%%%%%%%%%%%%%%%%%%%% REPRODUCE VALUES sl. 14
% Define ODE
dty = @(t,Y) t - 2*Y;                        % Derivative
y0  = 1;                                     % Initial condition

% Exact solution
Y = @(t) (1/4) * (2*t - 1 + 5*exp(-2*t));

% Numerical solutions
t0 = 0;
tf = 1;
Nt = 5;
[Tnum,Ynum] = odeSchemes(dty,[t0 tf],y0,Nt,'euler');

% Errors
err = Y(Tnum) - Ynum;
disp('Solution of the analytical problem slide 14 :')
disp([Tnum Ynum err])


%%%%%%%%%%%%%%%%%%%%%%%%% REPRODUCE CURVES sl. 15 & 19
% Analytic soution
t0   = 0;
tf   = 2;
Tref = linspace(t0,tf,1000);
Yref  = Y(Tref);

% Graphical representation
figure
plot(Tref,Yref)
hold on
xlabel('t')
ylabel('Y(t)')
grid on

% Numerical solutions
Nt  = 10 * (2.^(0:4))';
err = zeros(size(Nt));
stp = zeros(size(Nt));
for i = 1:length(err)
    % Euler scheme
    [Tnum,Ynum] = odeSchemes(dty,[t0 tf],y0,Nt(i),'euler');
    
    % Step
    stp(i) = Tnum(2) - Tnum(1);
    
    % Error
    err(i) = max(abs(Ynum - Y(Tnum)));
    
    % Graphical representation
    plot(Tnum,Ynum,'o--')
end

% Information
legend({'Analytic','Numeric'})
disp('Error in function of Nt')
disp([Nt stp err])
    
% Final error
figure
loglog(stp,stp.^1,'k-')
hold on
loglog(stp,err,'-o')
grid on
xlabel('Step h')
ylabel('max(error)')
title('Convergence comparison (log scale)')
legend({'f(x))= x^1','Euler'})


%%%%%%%%%%%%%%%%%%%%%%%%% MASS-SPRING sl. 24
% Physical parameters
k  = 1;
m  = 1;
c  = 0;
z0 = -1;

% Define EDO
dtX = @(t,X) X * [-c/m -k/m ; 1 0]';                        % Derivative
X0  = [0 z0];                                            % Initial condition

% Exact solution
z = @(t) z0 * cos(k/m*t);

% Analytic soution
t0   = 0;
tf   = 20;
Tref = linspace(t0,tf,1000);
Zref = z(Tref);

% Graphical representation
figure
plot(Tref,Zref)
hold on
xlabel('t')
ylabel('Y(t)')
grid on

% Numerical solution
[Tnum,Xnum] = odeSchemes(dtX,[t0 tf],X0,100,'euler');
plot(Tnum,Xnum(:,2),'o-r')
legend({'Analytic','Euler'})
title('Solution of a time harmonic oscillator')


%%%%%%%%%%%%%%%%%%%%%%%%% PREDATOR PREY sl. 25
% Physical parameters
alpha = 2;
beta  = 0.02;
delta = 0.0002;
gamma = 0.8;

% Define EDO
dtX = @(t,X) [alpha*X(1)-beta*X(1)*X(2) , delta*X(1)*X(2)-gamma*X(2)]                        % Derivative
X0  = [5000 100];                                            % Initial condition

% Numerical solution
t0 = 0;
tf = 30;
[Tnum,Xnum] = odeSchemes(dtX,[t0 tf],X0,100,'RK4');

% Graphical representation
figure
title('Solution of a predator prey system')
subplot(2,1,1)
plot(Tnum,Xnum(:,1))
xlabel('t')
ylabel('Prey population')
grid on
subplot(2,1,2)
plot(Tnum,Xnum(:,2))
xlabel('t')
ylabel('Predator population')
grid on


