source: Sophya/trunk/SophyaLib/NTools/difeq.h@ 1557

// This may look like C code, but it is really -*- C++ -*-
//  Simple differential equation solver 
//        E. Aubourg 1996-
// DAPNIA/SPP (Saclay) / CEA    LAL - IN2P3/CNRS  (Orsay)

#ifndef DIFEQ_H_SEEN
#define DIFEQ_H_SEEN

#include "machdefs.h"
#include "pexceptions.h"
#include "tvector.h"

namespace SOPHYA {

class GeneralFunction;

// <summary> fonction pour equadifs </summary>
// Une fonction utilisee pour les equations differentielles.
// On resoud de facon generale un systeme de n equations differentielles,
// donnant les derivees fpi de n fonctions fi.
class DiffEqFunction  {
public:
  // Constructeur. n = nombre de fonctions dans le systeme
  DiffEqFunction(int n)          : mNFunc(n), mNFuncApp(n) {}
  
  // Constructeur. n = nombre reel de fonctions dans le systeme,
  // m = nombre apparent (il y a dans certaines cas des fonctions a 
  // usage interne,  par exemple la fonction constante valant 1...)
  
  DiffEqFunction(int n, int m)   : mNFunc(n), mNFuncApp(m) {}
  
  // Destructeur
  virtual ~DiffEqFunction()      {}
  
  // Calcule les valeurs des derivees fpi a partir des valeurs des fonctions fi
  virtual void ComputeV(Vector& fpi, Vector const& fi)
       { Compute(fpi(0), fi(0)); }

  // Dans le cas ou on a une seule fonction, calcule la valeur de la derivee fp
  // a partir de la valeur de la fonction
  virtual void Compute(double& /*fp*/, double /*f*/)
       { ASSERT(mNFunc == 1); }
  
  // Nombre apparent de fonctions dans le systeme
  int NFunc()     {return mNFuncApp;}
  
  // Nombre reel de fonctions dans le systeme
  int NFuncReal() {return mNFunc;}
  
  // Pour ajuster vecteur de depart quand il y a des fonctions a usage
  // interne...
  virtual void AdjustStart(Vector& /*start*/, double /*tstart*/) 
    {}   
protected:
  // Nombre de fonctions dans le systeme
  int mNFunc;                   
  // Nombre apparent de fonctions                 
  int mNFuncApp;                               
};


// Cas y' = f(y)
typedef double(*DIFEQFCN1)(double);

// <summary> Cas y' = f(y) </summary>
// Cas y' = f(y), on fournit la fonction f, sous la forme
// double f(double), et ca construit la bonne DiffEqFunction
class DiffEqFcn1 : public DiffEqFunction {
public:
  // Constructeur, on fournit une fonction double->double
  // qui donne y' en fonction de y.
  DiffEqFcn1(DIFEQFCN1);

  // Implementation de Compute qui va utiliser la fonction fournie
  // au constructeur.
  virtual void Compute(double& fp, double f);
protected:
  DIFEQFCN1 mFcn;
};

// Cas y' = f(y,t)
typedef double(*DIFEQFCNT1)(double, double);

// <summary> y' = f(y,t) </summary>
// Cas y' = f(y,t), on fournit la fonction f, sous la forme
// double f(double, double), et ca construit la bonne DiffEqFunction
class DiffEqFcnT1 : public DiffEqFunction {
public:
  // Constructeur, on fournit une fonction (double, double)->double
  // qui donne y' en fonction de y et t.
  DiffEqFcnT1(DIFEQFCNT1);

  // Implementation de Compute qui va utiliser la fonction fournie
  // au constructeur.
  virtual void ComputeV(Vector& fpi, Vector const& fi);

  // Implementation de AdjustStart qui gere la fonction a usage interne.
  virtual void AdjustStart(Vector& start, double tstart);
protected:
  DIFEQFCNT1 mFcn;
};

// Cas y'' = f(y',y)
typedef double(*DIFEQFCN2)(double, double);

class DiffEqFcn2 : public DiffEqFunction {
public:
  DiffEqFcn2(DIFEQFCN2);

  virtual void ComputeV(Vector& fpi, Vector const& fi);
protected:
  DIFEQFCN2 mFcn;
};

// Cas y'' = f(y',y,t)
typedef double(*DIFEQFCNT2)(double, double, double);

// <summary> y'' = f(y',y,t) </summary>
// Cas y'' = f(y',y,t), on fournit la fonction f, sous la forme
// double f(double, double, double), et ca construit la bonne DiffEqFunction
class DiffEqFcnT2 : public DiffEqFunction {
public:
  // Constructeur, on fournit une fonction (double, double, double)->double
  // qui donne y'' en fonction de y', y et t.
  DiffEqFcnT2(DIFEQFCNT2);

  // Implementation de Compute qui va utiliser la fonction fournie
  // au constructeur.
  virtual void ComputeV(Vector& fpi, Vector const& fi);

  // Implementation de AdjustStart qui gere la fonction a usage interne.
  virtual void AdjustStart(Vector& start, double tstart);
protected:
  DIFEQFCNT2 mFcn;
};

// Cas y'' = f(y',y) avec des 3-vecteurs
typedef void(*DIFEQFCNV)(Vector&, Vector const&, Vector const&);

// <summary> y'' = f(y',y,t) </summary>
// Cas y'' = f(y',y,t), on fournit la fonction f, sous la forme
// double f(Vector), et ca construit la bonne DiffEqFunction
class DiffEqFcnV : public DiffEqFunction {
public:
  // Constructeur, on fournit une fonction (Vector)->double
  // qui donne y'' en fonction du vecteur (t, y, y')
  DiffEqFcnV(DIFEQFCNV);

  // Implementation de Compute qui va utiliser la fonction fournie
  // au constructeur.
  virtual void ComputeV(Vector& fpi, Vector const& fi);
protected:
  DIFEQFCNV mFcn;
  Vector tmp1, tmp2, tmp3;
};

// <summary> Classe abstraite de resolveur d'equadif </summary>
// Classe abstraite de resolveur d'equadif
// Beaucoup de fonctions renvoient l'objet pour pouvoir faire une
// notation chainee s.Step(...).Start(...).Solve(...)
class DiffEqSolver  {
public:
  // Constructeurs. L'equadif est donnee sous forme de DiffEqFunction.
  // On a prevu le cas particulier du premier degre directement
  // <group>
  DiffEqSolver();
  DiffEqSolver(DiffEqFunction*);
  DiffEqSolver(DIFEQFCN1);
  // </group>
  
  // Destructeur
  virtual ~DiffEqSolver();

  // Change la fonction. Notation chainee possible
  // <group>
  DiffEqSolver& Func(DiffEqFunction*);
  DiffEqSolver& Func(DIFEQFCN1);
  // </group>
  // Change le pas d'integration. Notation chainee possible
  DiffEqSolver& Step(double);
  
  // Change les conditions initiales. Notation chainee possible.
  // <group>
  DiffEqSolver& StartV(Vector const& yi, double t);
  // si NFunc == 1
  DiffEqSolver& Start1(double        yi, double t);
  DiffEqSolver& Start(double const* yi, double t);
  // </group>
  
  // Lance la resolution, avec ou sans conservation de n valeurs intermediaires
  // <group>
  virtual void SolveV(Vector& yf, double tf);
  // si NFunc == 1
  virtual void Solve1(double& yf, double tf);  
  virtual void Solve(double* yf, double tf);
  virtual void SolveArr(Matrix&  y, double* t, double tf, int n)=0;
  // si NFunc == 1
  virtual void SolveArr1(double*  y, double* t, double tf, int n);   
  virtual void SolveArr2(double** y, double* t, double tf, int n);
  // </group>
  
protected:
  DiffEqFunction* mFunc;
  bool mOwnFunc;

  Vector  mYStart;
  double  mXStart;
  double  mStep;
};

// <summary> Runge-Kutta ordre 4 </summary>
// Runge-Kutta ordre 4
class RK4DiffEq : public DiffEqSolver {
public:
  // Constructeurs. Voir <linkto class=DiffEqSolver>DiffEqSolver</linkto>
  // <group>
  RK4DiffEq();
  RK4DiffEq(DiffEqFunction*);
  RK4DiffEq(DIFEQFCN1);
  // </group>

  // Implementation de RK4
  virtual void SolveArr(Matrix& y, double* t, double tf, int n);

protected:
  // Un pas RK4
  void RKStep(Vector& newY, Vector const& y0, double dt);
  // Vecteurs utilises en interne, pour ne pas les reallouer.
  Vector k1, k2, k3, k4;
};


} // Fin du namespace

#endif