// This file is part of the Herezh++ application.
//
// The finite element software Herezh++ is dedicated to the field
// of mechanics for large transformations of solid structures.
// It is developed by Gérard Rio (APP: IDDN.FR.010.0106078.000.R.P.2006.035.20600)
// INSTITUT DE RECHERCHE DUPUY DE LÔME (IRDL) .
//
// Herezh++ is distributed under GPL 3 license ou ultérieure.
//
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
// AUTHOR : Gérard Rio
// E-MAIL : gerardrio56@free.fr
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty
// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see .
//
// For more information, please consult: .
/************************************************************************
* DATE: 12/02/2006 *
* $ *
* AUTEUR: G RIO (mailto:gerardrio56@free.fr) *
* $ *
* PROJET: Herezh++ *
* $ *
************************************************************************
* BUT: Algorithmes de base pour la résolution d'équations *
* différentielles. *
* $ *
* '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' *
* $ *
************************************************************************/
#ifndef ALGO_EDP_H
#define ALGO_EDP_H
#include "Vecteur.h"
#include "Mat_abstraite.h"
#include "Racine.h"
/** @defgroup Les_classes_algo
*
* BUT: Algorithmes de base pour la résolution d'équations
* différentielles, recherche de zéros, intégration.
*
*
* \author Gérard Rio
* \version 1.0
* \date 12/02/2006
* \brief Algorithmes de base pour la résolution d'équations différentielles, recherche de zéros, intégration.
*
*/
/// @addtogroup Les_classes_algo
/// @{
///
/// BUT: Algorithmes de base pour la résolution d'équations
/// différentielles.
class Algo_edp
{
public :
// CONSTRUCTEURS :
// constructeur par défaut
Algo_edp();
// constructeur de copie
Algo_edp(const Algo_edp& a);
// DESTRUCTEUR :
~Algo_edp();
// METHODES PUBLIQUES :
/// résolution par runge kutta imbriqué --> estimation d'erreur
/// il y a 3 calcul de la fonction dérivée
/// 2 et 3 ième ordre imbriquée
/// calcul de la solution à t+dt en fonction de la solution à t,
/// en entrée:
/// *Ptder_fonc : pointeur de la fonction de calcul de la dérivées
/// en entrée: t et f : temps et la valeur de la fonction a t
/// en sortie: df : dérivée de la fonction a t
/// erreur : si diff de 0, indique qu'il y a eu une erreur
/// *Ptverif_fonc : pointeur de fonction, pour tester val_finale après la prédiction explicite finale
///
/// val_initiale : valeur des fonctions en t0
/// der_initiale : valeur de la dérivée pour t0
/// t0 : temps initiale
/// deltat : incrément de temps demandé
/// en sortie :
/// val_finale : valeur des fonctions à tdt
/// estime_erreur : estimation d'erreur pour chaque fonction
/// --> si estime_erreur(i) est >= à ConstMath::tresgrand, la valeur finale = la valeur initiale (c-a-d pas de calcul valide)
template void Runge_Kutta_step23(T& instance
,Vecteur& (T::*Ptder_fonc) (const double & t,const Vecteur& f,Vecteur& df,int& erreur)
,void (T::*Ptverif_fonc) (const double & t,const Vecteur& f,int& erreur_final)
,const Vecteur& val_initiale,const Vecteur& der_initiale
,const double& t0,const double& deltat
,Vecteur& val_finale,Vecteur& estime_erreur) ;
/// résolution par runge kutta imbriqué --> estimation d'erreur
/// il y a 4 calcul de la fonction dérivée
/// 3 et 4 ième ordre imbriquée (Fehlberg)
/// calcul de la solution à t+dt en fonction de la solution à t,
/// en entrée:
/// *Ptder_fonc : pointeur de la fonction de calcul de la dérivées
/// en entrée: t et f : temps et la valeur de la fonction a t
/// en sortie: df : dérivée de la fonction a t
/// erreur : si diff de 0, indique qu'il y a eu une erreur
/// *Ptverif_fonc : pointeur de fonction, pour tester val_finale après la prédiction explicite finale
///
/// val_initiale : valeur des fonctions en t0
/// der_initiale : valeur de la dérivée pour t0
/// t0 : temps initiale
/// deltat : incrément de temps demandé
/// en sortie :
/// val_finale : valeur des fonctions à tdt
/// estime_erreur : estimation d'erreur pour chaque fonction
/// --> si estime_erreur(i) est >= à ConstMath::tresgrand, la valeur finale = la valeur initiale (c-a-d pas de calcul valide)
template void Runge_Kutta_step34(T& instance
,Vecteur& (T::*Ptder_fonc) (const double & t,const Vecteur& f,Vecteur& df,int& erreur)
,void (T::*Ptverif_fonc) (const double & t,const Vecteur& f,int& erreur_final)
,const Vecteur& val_initiale,const Vecteur& der_initiale
,const double& t0,const double& deltat
,Vecteur& val_finale,Vecteur& estime_erreur) ;
/// résolution par runge kutta imbriqué --> estimation d'erreur
/// il y a 5 calcul de la fonction dérivée
/// l'algorithme s'appuis sur la présentation de numerical recipes
/// fifth-order Cash-Karp Runge-Kutta
/// calcul de la solution à t+dt en fonction de la solution à t,
/// en entrée:
/// *Ptder_fonc : pointeur de la fonction de calcul de la dérivées
/// en entrée: t et f : temps et la valeur de la fonction a t
/// en sortie: df : dérivée de la fonction a t
/// erreur : si diff de 0, indique qu'il y a eu une erreur
/// *Ptverif_fonc : pointeur de fonction, pour tester val_finale après la prédiction explicite finale
///
/// val_initiale : valeur des fonctions en t0
/// der_initiale : valeur de la dérivée pour t0
/// t0 : temps initiale
/// deltat : incrément de temps demandé
/// en sortie :
/// val_finale : valeur des fonctions à tdt
/// estime_erreur : estimation d'erreur pour chaque fonction
/// --> si estime_erreur(i) est >= à ConstMath::tresgrand, la valeur finale = la valeur initiale (c-a-d pas de calcul valide)
template void Runge_Kutta_step45(T& instance
,Vecteur& (T::*Ptder_fonc) (const double & t,const Vecteur& f,Vecteur& df,int& erreur)
,void (T::*Ptverif_fonc) (const double & t,const Vecteur& f,int& erreur_final)
,const Vecteur& val_initiale,const Vecteur& der_initiale
,const double& t0,const double& deltat
,Vecteur& val_finale,Vecteur& estime_erreur) ;
/// un programme de pilotage de l'intégration d'un pas (on s'inspire du programme libre rkf45)
/// on distingue une erreur globale et une erreur relative
/// l'erreur réelle utilisée sur chaque step est: erreurRelative * |f| + erreurAbsolue
/// ==== en entrée: =====
/// cas_kutta : donne le type de routine kutta imbriqué a employer:
/// =3 pour kutta23, =4 pour kutta34, =5 pour kutta45
/// tdebut et tfin : temps de début et de fin du calcul, demandé
/// val_initiale : valeur initiale de la fonction (donc à tdebut)
/// der_initiale : dérivée initiale de la fonction (donc à tdebut)
/// erreurAbsolue : erreur absolue demandée
/// erreurRelative : erreur relative (à la fonction) demandée
/// Ptder_fonc : pointeur de fonction
/// en entrée: t et f : temps et la valeur de la fonction a t
/// en sortie: df : dérivée de la fonction a t
/// en entrée: t et f : temps et la valeur de la fonction a t
/// en sortie: df : dérivée de la fonction a t
/// erreur : si diff de 0, indique qu'il y a eu une erreur
/// Ptverif_fonc : pointeur de fonction, pour tester l'intégrité du résultat
/// cela signifie que val_finale est testée avec cette fonction, systématiquement
/// - à la fin de chaque prédiction des kuttas imbriqués
/// - au retour du pilotage
/// ---- en sortie: -----
/// dernierTemps : temps final utilisé = tfin si calcul ok, sinon = le dernier temps calculé
/// val_finale : valeur de la fonction à "dernierTemps"
/// der_finale : dérivée finale
/// dernierdeltat : le dernier deltat utilisé
/// nombreAppelF : nombre d'appel de la fonction réellement utilisé dans le programme
/// erreurRelative : si retour = 3 ==> erreur relative (à la fonction) minimum possible,
/// nb_step : nombre de step utilisé pour le calcul (non compté les steps avec trop d'erreur)
/// erreur_maxi_global : une approximation de l'erreur maxi global cumulant les erreurs sur tous les steps
/// === en retour un entier qui donne le résultat du calcul : ====
/// =2 : tout est ok : seul cas où toutes les valeurs de retour sont valides
/// =3 : erreur: la précision relative demandée est trop petite, en retour la précision mini possible
/// =4 : erreur: on a dépassé le nombre maxi d'appel fixé, d'où a peu près à (nombreAppelF/6) step de calcul.
/// =6 : erreur: l'intégration pas possible, due aux précisions demandées, on doit augmenter ces précisions
/// souvent du à une solution qui varie très rapidemment --> problème potentiel
/// =8 : erreur: cas_kutta ne correspond pas à une routine de base existante
/// = 0 : erreur inconnue
template
int Pilotage_kutta(int cas_kutta, T& instance
,Vecteur& (T::*Ptder_fonc) (const double & t,const Vecteur& f,Vecteur& df,int& erreur)
,void (T::*Ptverif_fonc) (const double & t,const Vecteur& f,int& erreur_final)
,const Vecteur& val_initiale,const Vecteur& der_initiale
,const double& tdebut,const double& tfin
,const double& erreurAbsolue, double& erreurRelative
,Vecteur& val_finale,Vecteur& der_finale
,double& dernierTemps, double& dernierdeltat,int& nombreAppelF,int& nb_step
,double& erreur_maxi_global);
/// explication erreur
/// affiche à l'écran l'explication de l'erreur de Pilotage_kutta
void Affiche_Explication_erreur_RG(int type_erreur);
/// --méthodes pour modifier les différents paramètres
/// nombre maxi d'appel de fonction permis
/// dans le cas du runge 4-5 -> 6 appels par step (à la louche)
void Modif_nbMaxiAppel (int nb) {nbMaxiAppel = nb;};
/// initialisation de tous les paramètres à leurs valeurs par défaut
void Init_param_val_defaut();
/// modifie le niveau d'affichage par exemple pour le débug
void Change_niveau_affichage(int niveau) {permet_affichage=niveau;};
private :
// VARIABLES PROTEGEES :
// ----- controle de la sortie des informations
int permet_affichage; // pour permettre un affichage spécifique dans les méthodes, pour les erreurs et warning
int nbMaxiAppel; // nombre maxi d'Appel autorisé de la fonction
// ---- variables intermédiaires utilisées par la méthode Runge_Kutta_step23
static double AB1,AB2,AB3,A2,aa2,aa3
,bb21,bb31,bb32
,dA1,dA2,dA3;
// ---- variables intermédiaires utilisées par la méthode Runge_Kutta_step34
static double A_B1,A_B3,A_B4,A_B5,A_1,A_3,A_4
,a_a2,a_a3,a_a4,a_a5
,b_b21,b_b31,b_b32,b_b41,b_b42,b_b43,b_b51,b_b53,b_b54
,d_A1,d_A3,d_A4,d_A5;
// --- variables intermédiaires utilisées par la méthode Runge_Kutta_step45
Vecteur k2,k3,k4,k5,k6;
// Vecteur f_inter;
static double a2,a3,a4,a5,a6 // les ai de Cash-karp paramètres
,b21,b31,b32,b41,b42,b43 //
,b51,b52,b53,b54,b61,b62,b63,b64,b65 // les bij
,c1,c3,c4,c6; // les ci
// c2 =0, c5=0
static double ce1,ce3,ce4,ce5,ce6; // ce2=0; : les c-c_étoile_i
// --- variables intermédiaires utilisées par la méthode Pilotage_kutta
Vecteur val_inter; // valeur de la fonction, intermédiaire
Vecteur der_inter; // dérivée de la fonction, intermédiaire
Vecteur estime_erreur; // estimation d'erreur calculée
// puis des sauvegardes pour pouvoir revenir en arrière
Vecteur val_prec_inter; // sauve avant appel: valeur de la fonction, intermédiaire
Vecteur der_prec_inter; // sauve avant appel: dérivée de la fonction, intermédiaire
Vecteur estime_prec_erreur; // sauve avant appel: estimation d'erreur calculée
double coef; // coeff de sécurité préconnisé pour les variations de temps
static double remin; // pour calculer l'erreur mini acceptable dans le runge
// METHODES PROTEGEES :
// impression du temps
void timestamp ();
// -- méthode test ---
// on défini une fonction a intégrer
Vecteur& FoncDeriv(const double & t,const Vecteur& f,Vecteur& df,int& erreur);
// on défini une fonction de vérification d'intégrité
void FoncVerif_fonc(const double & t,const Vecteur& f,int& erreur);
int nb_cas_test; // nb pour choisir entre différents cas test
// la méthode qui fait les tests
void Test_template();
};
// pour faire de l'inline: nécessaire avec les templates
// on n'inclut que les méthodes templates
#include "Algo_edp_2.cc"
#define Algo_edp_deja_inclus
/// @} // end of group
#endif