% Analyzing the accuracy of various schemes 
% for numerical solution of ODEs

function ode_02
  
echo off, clc
clear all
close all

% Code options;
STABILITY = 1;
ACCURACY  = 0;
CONSERVATION = 0;

% Parameters of the equation;
global Lr;  Lr = -1.0; 
global Li;  Li = 1.0 * j;
global y0;  y0 = 1.0;

% analyzing stability
if ( STABILITY )
% Length of the integration interval;
  T = 100.0;
  
  Euler_Explicit( 0.05, T, 1 ); 
  Euler_Explicit( 2.0, T, 2 ); 

%  Crank_Nicholson( 1.0, T, 1 ); 
%  Crank_Nicholson( 2.05, T, 2 ); 

%  Leapfrog( 0.05, T, 1 );
%  Leapfrog( 1.05, T, 2 ); 

end;



% analyzing accuracy
if ( ACCURACY ) 
% Length of the integration interval;
  T = 10.0;
  N = 10;
  Nmax = 1e4;
  h = T / N;
  
  i = 1;
  while ( N <= Nmax );  
    h = T / N;
    [E1r E1i] = Euler_Explicit( h, T, 0 );
    [E2r E2i] = Euler_Implicit( h, T, 0 );
    [E3r E3i] = Crank_Nicholson( h, T, 0 );
    [E4r E4i] = Leapfrog( h, T, 0 );
    [E5r E5i] = RK4( h, T, 0 );
  
    Err1(i,1) = h;
    Err2(i,1) = h;
    Err1(i,2) = E1r;
    Err2(i,2) = E1i;
    Err1(i,3) = E2r;
    Err2(i,3) = E2i;
    Err1(i,4) = E3r;
    Err2(i,4) = E3i;
    Err1(i,5) = E4r;
    Err2(i,5) = E4i;
    Err1(i,6) = E5r;
    Err2(i,6) = E5i;

    i = i + 1;
    N = N * 2;
  end;
  
  figure(1); loglog(Err1(:,1),Err1(:,2),'-x',Err1(:,1),Err1(:,3),'-+',Err1(:,1),Err1(:,4),'-o',Err1(:,1),Err1(:,5),'-^',Err1(:,1),Err1(:,6),'-d');
  legend('Explicit Euler','Implicit Euler','Crank-Nicholson','Leapfrog','RK4'); 
  title('Real');
  xlabel('h'); ylabel('error');

  figure(2); loglog(Err2(:,1),Err2(:,2),'-x',Err2(:,1),Err2(:,3),'-+',Err2(:,1),Err2(:,4),'-o',Err2(:,1),Err2(:,5),'-^',Err2(:,1),Err2(:,6),'-d');
  legend('Explicit Euler','Implicit Euler','Crank-Nicholson','Leapfrog','RK4'); 
  title('Imag');
  xlabel('h'); ylabel('error');
end;


% Integration of the model problem using the Euler Explicit scheme;
function [eT1, eT2] = Euler_Explicit( h, T, fn )
  
global Lr; global Li; global y0;

i = 1;
t = i * h;
Yrn = y0;
Yin = y0;
while ( t <= T )
  t = i * h;
  Yrn1 = Yrn + h * (Lr * Yrn);
  Yin1 = Yin + h * (Li * Yin);
  Yrn = Yrn1;
  Yin = Yin1;
  Y(i,1) = t;
  Y(i,2) = Yrn1;
  Y(i,3) = y0 * exp( Lr * t );
  Y(i,4) = imag(Yin1);
  Y(i,5) = imag(y0 * exp( Li * t ));  
  i = i + 1;
end;

if ( fn > 0 )
  figure(fn); 
  subplot(2,1,1);
  plot(Y(:,1),Y(:,2),'-b',Y(:,1),Y(:,3),'-.b'); title('Real');
  subplot(2,1,2);
  plot(Y(:,1),Y(:,4),'-r',Y(:,1),Y(:,5),'-.r'); title('Imag');
end;

eT1 = abs( Y(i-1,2) - Y(i-1,3) ) / abs(Y(i-1,3));
eT2 = abs( Y(i-1,4) - Y(i-1,5) ) / abs(Y(i-1,5));




% Integration of the model problem using the Euler Implicit scheme;
function [eT1, eT2] = Euler_Implicit( h, T, fn )
  
global Lr; global Li; global y0;

i = 1;
t = i * h;
Yrn = y0;
Yin = y0;
while ( t <= T )
  t = i * h;
  Yrn1 = Yrn / ( 1.0 - h * Lr );
  Yin1 = Yin / ( 1.0 - h * Li );
  Yrn = Yrn1;
  Yin = Yin1;
  Y(i,1) = t;
  Y(i,2) = Yrn1;
  Y(i,3) = y0 * exp( Lr * t );
  Y(i,4) = imag(Yin1);
  Y(i,5) = imag(y0 * exp( Li * t ));  
  i = i + 1;
end;

if ( fn > 0 )
  figure(fn); 
  subplot(2,1,1);
  plot(Y(:,1),Y(:,2),'-b',Y(:,1),Y(:,3),'-.b'); title('Real');
  subplot(2,1,2);
  plot(Y(:,1),Y(:,4),'-r',Y(:,1),Y(:,5),'-.r'); title('Imag');
end;

eT1 = abs( Y(i-1,2) - Y(i-1,3) ) / abs(Y(i-1,3));
eT2 = abs( Y(i-1,4) - Y(i-1,5) ) / abs(Y(i-1,5));


% Integration of the model problem using the Crank-Nicholson scheme;
function [eT1, eT2] = Crank_Nicholson( h, T, fn )
  
global Lr; global Li; global y0;

i = 1;
t = i * h;
Yrn = y0;
Yin = y0;
while ( t <= T )
  t = i * h;
  Yrn1 = Yrn * ( 1.0 + 0.5 * h * Lr ) / ( 1.0 - 0.5 * h * Lr );
  Yin1 = Yin * ( 1.0 + 0.5 * h * Li ) / ( 1.0 - 0.5 * h * Li );
  Yrn = Yrn1;
  Yin = Yin1;
  Y(i,1) = t;
  Y(i,2) = Yrn1;
  Y(i,3) = y0 * exp( Lr * t );
  Y(i,4) = imag(Yin1);
  Y(i,5) = imag(y0 * exp( Li * t ));  
  i = i + 1;
end;

if ( fn > 0 )
  figure(fn); 
  subplot(2,1,1);
  plot(Y(:,1),Y(:,2),'-b',Y(:,1),Y(:,3),'-.b'); title('Real');
  subplot(2,1,2);
  plot(Y(:,1),Y(:,4),'-r',Y(:,1),Y(:,5),'-.r'); title('Imag');
end;

eT1 = abs( Y(i-1,2) - Y(i-1,3) ) / abs(Y(i-1,3));
eT2 = abs( Y(i-1,4) - Y(i-1,5) ) / abs(Y(i-1,5));


% Integration of the model problem using the Leapfrog scheme;
function [eT1, eT2] = Leapfrog( h, T, fn )
  
global Lr; global Li; global y0;

i = 2;
t = i * h;
Yrna = y0;
Yina = y0;
Yrnb = y0 * exp( Lr * h );
Yinb = y0 * exp( Li * h );

while ( t <= T )
  t = i * h;
  if ( mod(i,2) == 0 )
    Yrn1a = Yrna + 2.0 * h * Lr * Yrnb;
    Yin1a = Yina + 2.0 * h * Li * Yinb;
    Yrna = Yrn1a;
    Yina = Yin1a;
    Y(i,2) = Yrn1a;
    Y(i,4) = imag(Yin1a);
  else
    Yrn1b = Yrnb + 2.0 * h * Lr * Yrna;
    Yin1b = Yinb + 2.0 * h * Li * Yina;
    Yrnb = Yrn1b;
    Yinb = Yin1b;
    Y(i,2) = Yrn1b;
    Y(i,4) = imag(Yin1b);    
  end;
  Y(i,1) = t;
  Y(i,3) = y0 * exp( Lr * t );
  Y(i,5) = imag(y0 * exp( Li * t ));  
  i = i + 1;
end;

if ( fn > 0 )
  figure(fn); 
  subplot(2,1,1);
  plot(Y(:,1),Y(:,2),'-b',Y(:,1),Y(:,3),'-.b'); title('Real');
  subplot(2,1,2);
  plot(Y(:,1),Y(:,4),'-r',Y(:,1),Y(:,5),'-.r'); title('Imag');
end;

eT1 = abs( Y(i-1,2) - Y(i-1,3) ) / abs(Y(i-1,3));
eT2 = abs( Y(i-1,4) - Y(i-1,5) ) / abs(Y(i-1,5));



% Integration of the model problem using the RK4 scheme;
function [eT1, eT2] = RK4( h, T, fn )
  
global Lr; global Li; global y0;

i = 1;
t = i * h;
Yrn = y0;
Yin = y0;
while ( t <= T )
  t = i * h;
  k1r = Lr * Yrn;
  k1i = Li * Yin;
  
  k2r = Lr * (Yrn + 0.5 * h * k1r);
  k2i = Li * (Yin + 0.5 * h * k1i);
  
  k3r = Lr * (Yrn + 0.5 * h * k2r);
  k3i = Li * (Yin + 0.5 * h * k2i);
  
  k4r = Lr * (Yrn + h * k3r);
  k4i = Li * (Yin + h * k3i);
  
  Yrn1 = Yrn + h * ( k1r/6.0 + (k2r + k3r)/3.0 + k4r/6.0 );
  Yin1 = Yin + h * ( k1i/6.0 + (k2i + k3i)/3.0 + k4i/6.0 );
  Yrn = Yrn1;
  Yin = Yin1;
  Y(i,1) = t;
  Y(i,2) = Yrn1;
  Y(i,3) = y0 * exp( Lr * t );
  Y(i,4) = imag(Yin1);
  Y(i,5) = imag(y0 * exp( Li * t ));  
  i = i + 1;
end;

if ( fn > 0 )
  figure(fn); 
  subplot(2,1,1);
  plot(Y(:,1),Y(:,2),'-b',Y(:,1),Y(:,3),'-.b'); title('Real');
  subplot(2,1,2);
  plot(Y(:,1),Y(:,4),'-r',Y(:,1),Y(:,5),'-.r'); title('Imag');
end;

eT1 = abs( Y(i-1,2) - Y(i-1,3) ) / abs(Y(i-1,3));
eT2 = abs( Y(i-1,4) - Y(i-1,5) ) / abs(Y(i-1,5));
