%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% COURS: STATGE DE RENTREE
%
% @copyright: J.-B.A.K. <jean-baptiste.apoung@math.u-psud.fr>
% ANNEE : 2015-2016
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [t,x,xe,res] = AdaptiveOneStep()
  
  t0 = 0;
  x0 = 1;
  tf = pi/3;
  eps=1e-4;
  
  %decommenter pour executer
  
  %theme1();
  %[t,x,xe] = monTest(t0,tf,x0);
  %plot(t,x,'*', t,xe,'-');
  [t,x,xe] = monTestNaif(t0,tf,x0,eps);figure(1)
  plot(t,x,'*', t,xe,'-');
  [t,x,xe] = monTestEvolue(t0,tf,x0,10 * eps); figure(2)
  plot(t,x,'*', t,xe,'-');
  [t,x,xe] = monTestEvolueApprochee(t0,tf,x0,eps/100); figure(3)
  plot(t,x,'*', t,xe,'-');
  res=1;
end
%NERIEN TOUCHER APRES CETTE LIGNE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Theme 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function []=theme1()
  t0=0;
  tf = 6.6632868593231301896996820305;
  A  = 2.00861986087484313650940188;
  x0=[A;0];
  
  mu=[10^-1, 10^0, 10^1, 10^2, 10^3];
  nbode45=[];
  nbode23=[];
  for i =1:1:3
    f = @(t,x) VanderPol(t,x,mu(i));
    [t,x]=ode45(f, [t0,tf],x0);
    nbode45=[nbode45, length(t)];
    [t,x]=ode23(f, [t0,tf],x0);
    nbode23=[nbode23, length(t)];
  end
  figure(1)
  plot(x(:,1),x(:,2),'+');
  figure(2)
  plot(mu,nbode45, '-+r', mu,nbode23,'-*b');
  legend('ode45','odeq23');
end
%%------------------------------------------------------------------------------------
function [f, ft, fy] = VanderPol(t,x,mu)
  % f  est  la valeur de la fonction second-membre f  au point (t,y)
  % ft est la dériée partielle parrapport a t de la de f au point (t,y) (i.e df/dt )
  % fy est la dériée partielle parrapport a t de la de f au point (t,y) (i.e df/dy )
  f=zeros(2,1);
  f(1) = x(2);
  f(2) = mu *(1 -  x(1)*x(1))*x(2) - x(1) ;
  fy  = [0, 1; mu*2*x(1)*x(2)-1, mu*(1 -  x(1)*x(1))];
  ft = [0, 0;];
end
%%-------------------------------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Theme 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
% Exercice 2
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function [t,x,xe] = monTest(t0,tf,x0)
  t=t0;
  x = x0;
  format='%9.9e  %9.9e  %9.9e\n';
    ordre=2;
    h = 1/1000;
    xe = solExacte(t(end));
    err = 0;
    while (t(end) < tf)
      h = min(h, tf - t(end));
      [xx,errr] = monUnPasErreur(t(end),x(end),err,h);
      x=[x,xx]; err = errr;
      t=[t,t(end)+h];
      xee = solExacte(t(end));
      error = abs(xee - xx);
      xe=[xe,xee];
      printf(format, t(end), error, abs(errr*h*h));
    end
  end
  %end-monTestSimple
  %--------------------------------------------------------
  function [xnn,enn] = monUnPasErreur(tn,xn,en,h)
    %% xnn correspond  a xn+1
    %% enn correspond  a en+1
    %[phi, tau, fy] = monRK2(tn, xn, h);
    [f, ft, fy] = monEDO (tn,xn);
    p1 = f;
    p2 = monEDO(tn + 2 * h/3, xn + 2 * h * p1/3);
    phi = (p1 + 3 * p2)/4 ;
    tau = fy * (ft + fy * f) / 6;
    xnn = xn + h * phi;
    enn = en + h * ( fy * en + tau );
  end
  %end-monUnPasNaif
  %------------------------------------------------------------
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Theme 3
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  % Exercice 1
  %!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  %---------------------------------------------------------------------
  function [xnn,enn] = monUnPasNaif(tn,xn,en,h)
    %% xnn correspond  a xn+1
    %% enn correspond  a en+1
    %[phi, tau, fy] = monRK2(tn, xn, h);
    [f, ft, fy] = monEDO (tn,xn);
    p1 = f;
    p2 = monEDO(tn + 2 * h/3, xn + 2 * h * p1/3);
    phi = (p1 + 3 * p2)/4 ;
    tau = fy * (ft + fy * f) / 6;
    xnn = xn + h * phi;
    enn = en + h * ( fy * en + tau );
  end %monUnPasNaif
  %---------------------------------------------------------------------
  function [t,x,xe] = monTestNaif(t0,tf,x0,eps)
    t = t0;
    x = x0;
    xe = solExacte(t(end));
    ordre = 2;
    h = pasInitial(t0,x0,@monEDO,eps,ordre)
    t_courant = t0;
    x_courant = x;
    e_courant = 0;
    hmin=(tf-t0)/1e12;
    hmax = (tf-t0)/2.5;
    while (t_courant < tf)
      h = max(hmin, h);
      h = min(hmax,h);
      h = min(h, tf - t_courant);
      [xx, ee] = monUnPasNaif(t_courant,x_courant,e_courant,h) ;
      ratio = abs(ee)*(h^ordre)/(eps);
      if(ratio < 1)
        x_courant = xx;
        e_courant = ee;
        t_courant = t_courant + h;
        x=[x,x_courant];
        t=[t,t_courant];
        xe = [xe,solExacte(t(end))];
      end
      h =  max(0.2 * h, min(5.0 * h, (0.8 * eps/abs(ee))^(1/ordre)));
    end
  end
  %!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  % Exercice 2
  %!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  %-------------------------------------------
  function [xnn,tau] = monUnPasEvolue(tn,xn,h)
    % xnn correspond  a xn+1
    %[phi, tau, fy] = monRK2(tn, xn, h);
    [f, ft, fy] = monEDO (tn,xn);
    p1 = f;
    p2 = monEDO(tn + 2 * h/3, xn + 2 * h * p1/3);
    phi = (p1 + 3 * p2)/4;
    tau = fy * (ft + fy * f) / 6;
    xnn = xn + h * phi;
  end %monUnPasEVolue
  %-----------------------------------------
  function [t,x,xe] = monTestEvolue(t0,tf,x0,eps)
    t = t0;
    x = x0;
    xe = solExacte(t(end));
    ordre = 2;
    h = pasInitial(t0,x0,@monEDO,eps,ordre)
    t_courant = t0;
    x_courant = x;
    hmin=(tf-t0)/1e12;
    hmax = (tf-t0)/2.5;
    while (t_courant < tf)
      h = max(hmin, h);
      h = min(hmax,h);
      h = min(h, tf - t_courant);
      [xx, tau] = monUnPasEvolue(t_courant,x_courant,h) ;
      ratio = abs(tau)*(h^ordre)/(eps);
      if(ratio < 1)
        x_courant = xx;
        t_courant = t_courant + h;
        x=[x,x_courant];
        t=[t,t_courant];
        xe = [xe,solExacte(t(end))];
      end
      h =  max(0.2 * h, min(5.0 * h, (0.8 * eps/abs(tau))^(1/ordre)));
    end
  end %monTestEvolue
  
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Theme 4
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  function [t,x,xe] = monTestEvolueApprochee(t0,tf,x0,eps)
    t = t0;
    x = x0;
    xe = solExacte(t(end));
    ordre = 4;
    h = pasInitial(t0,x0,@monEDO,eps,ordre)
    t_courant = t0;
    x_courant = x;
    hmin=(tf-t0)/1e12;
    hmax = (tf-t0)/2.5;
    while (t_courant < tf)
      h = max(hmin, h);
      h = min(hmax,h);
      h = min(h, tf - t_courant);
      [xx, tau, ordre] = monUnPasDopri5(t_courant,x_courant,h,@monEDO) ;
      %  [xx, tau, ordre] =monUnPasFehlberg45(t_courant,x_courant,h,@monEDO) ;
      ratio = abs(tau)*(h^ordre)/(eps);
      if(ratio < 1)
        x_courant = xx;
        t_courant = t_courant + h;
        x=[x,x_courant];
        t=[t,t_courant];
        xe = [xe,solExacte(t(end))];
      end
      %    if (ratio < 1e-6) ratio = 1e-6; end
      %    coeff = max(0.2, min(5.0, ratio^(-1/ordre)));
      %    h = 0.8 * coeff * h;
      h =  max(0.2 * h, min(5.0 * h, (0.8 * eps/abs(tau))^(1/ordre)));
      
    end
  end
  %------------------------------------------------------------------------------------------
  function [xnn, tau, ordre] = monUnPasDopri5(t,x,h,f)
    % xnn -> solution calculee
    % tau -> approximation de la fonction principale d'erreur
    % p   -> ordre du schéma
    % f   est le second membre de l'equation differentielle a resoudre
    %%REMARQUE: pour plus d'efficacite on aurai pu retourner plutot ynn a la
    %%          place de tau, car on a la relation  tau = (ynn-xnn)/h^(ordre+1)
    k1 = f(t,x);
    k2 = f(t+h/4, x + h*k1/4);
    k3 = f(t+3*h/8, x+h*(3*k1/32 + 9*k2/32));
    k4 = f(t+h*12/13, x+ h*(1932*k1/2197 - 7200*k2/2197+ 7296*k3/2197));
    k5 = f(t+h, x+h*(439*k1/216 -8*k2 +3680*k3/513 -845*k4/4104));
    k6 = f(t+h/2, x+h*(-8*k1/27 + 2*k2 -3544*k3/2565 + 1859*k4/4104 -11*k5/40));
    xnn  = x + h*(25*k1/216 + 1408*k3/2565 + 2197*k4/4104 - k5/5);
    ynn  = x + h*(16*k1/135 + 6656*k3/12825 +  28561*k4/56430 - 9*k5/50 + 2*k6/55);
    tau  = (ynn - xnn)/(h^5);
    ordre = 4;
  end
  %------------------------------------------------------------------------------------------
  function [xnn,tau,ordre] = monUnPasFehlberg45(t,x,h,f)
    k1 = f(t,x);
    k2 = f(t+h/4, x + h*k1/4);
    k3 = f(t+3*h/8, x+h*(3*k1/32 + 9*k2/32));
    k4 = f(t+h*12/13, x+ h*(1932*k1/2197 - 7200*k2/2197+ 7296*k3/2197));
    k5 = f(t+h, x+h*(439*k1/216 -8*k2 +3680*k3/513 -845*k4/4104));
    k6 = f(t+h/2, x+h*(-8*k1/27 + 2*k2 -3544*k3/2565 + 1859*k4/4104 -11*k5/40));
    xnn = x + h*(25*k1/216 + 1408*k3/2565 + 2197*k4/4104 - k5/5);
    ynn = x + h*(16*k1/135 + 6656*k3/12825 +  28561*k4/56430 - 9*k5/50 + 2*k6/55);
    tau  = (ynn - xnn)/(h^5);
    ordre = 4;
  end
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%
  %%  UTILITAIRES
  %%
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  function [h] =pasInitial(t0,x0,f,tol,ordre)
    % fontion qui retourne le pas initial
    % pour resoudre une edo par un schema de pas  addaptatif
    %
    % t0  -> instant initial
    % x0  -> solution initial
    % f   -> second membre de l'EDO
    % tol -> tolerance à l'erreur
    % ordre -> l'ordre du schema considéré
    atol = tol;
    rtol = tol;
    n=length(x0);
    f1 = f(t0,x0);
    tol = atol + rtol * sqrt( dot(x0,x0)/n);
    d0 = sqrt( dot(x0,x0)/n)/tol;
    d1 = sqrt( dot(f1,f1)/n)/tol;
    if( (d0<1e-5) || (d1 < 1e-5))
      h0 = 1e-6;
    else
      h0 = 0.01*(d0/d1);
    end
    x1 = x0 + h0*f(t0,x0);
    f2 = f(t0+h0,x1);
    vd2=f2 -f1;
    d2 = sqrt( dot(vd2,vd2)/n)/(tol*h0);
    r=max(d0,d1);
    if(r <1e-15)
      h1 = max(1e-6,h0*1e-3);
    else
      h1 = (0.01/max(d1,d2))^(1.0/(ordre));
    end
    h = min(100*h0, h1);
  end
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  function [f, ft, fy] = monEDO (t,y)
    % f  est  la valeur de la fonction second-membre f  au point (t,y)
    % ft est la dériée partielle parrapport a t de la de f au point (t,y) (i.e df/dt )
    % fy est la dériée partielle parrapport a t de la de f au point (t,y) (i.e df/dy )
    f = (cos(t) - 2 * t * tan(t*t))*y;
    ft = -y*(sin(t) + 2 * tan(t*t) + 4 * t*t*(1+ tan(t*t)^2) );
    fy = cos(t) - 2*t*tan(t*t);
  end
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  function x = solExacte(t)
    x = exp(sin(t))*cos(t^2);
  end
