620 lines
32 KiB
C++
620 lines
32 KiB
C++
|
// 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) <https://www.irdl.fr/>.
|
||
|
|
//
|
||
|
|
// Herezh++ is distributed under GPL 3 license ou ultérieure.
|
||
|
|
//
|
||
|
|
// Copyright (C) 1997-2021 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 <https://www.gnu.org/licenses/>.
|
||
|
|
//
|
||
|
|
// For more information, please consult: <https://herezh.irdl.fr/>.
|
||
|
|
|
||
|
|
#ifndef Algo_edp_deja_inclus
|
||
|
|
#include "Algo_edp.h"
|
||
|
|
#endif
|
||
|
|
|
||
|
|
// --- include pour rkf45 ---
|
||
|
|
# include <cstdlib>
|
||
|
|
# include <iostream>
|
||
|
|
# include <iomanip>
|
||
|
|
# include <cmath>
|
||
|
|
# include <ctime>
|
||
|
|
// --- fin include pour rkf45 ---
|
||
|
|
|
||
|
|
// contient que les méthodes nécessitant les templates, qui doivent donc être
|
||
|
|
//incluses dans le .h
|
||
|
|
|
||
|
|
#ifndef Algo_edp_deja_inclus
|
||
|
|
#include "MathUtil.h"
|
||
|
|
#include "ConstMath.h"
|
||
|
|
static double remin = 1.0E-12;
|
||
|
|
|
||
|
|
|
||
|
|
// ----- cas d'un système d'edp du premier ordre---- ,
|
||
|
|
// 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
|
||
|
|
template <class T>
|
||
|
|
void Algo_edp::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écup de la dimension
|
||
|
|
int taille = val_initiale.Taille();
|
||
|
|
int erreur=0; // init de l'indicateur d'erreur dans le calcul de la fonction
|
||
|
|
// vérif de la dimmension
|
||
|
|
if (taille != k2.Taille())
|
||
|
|
{ k2.Change_taille(taille);k3.Change_taille(taille);
|
||
|
|
// val_finale.Change_taille(taille);
|
||
|
|
};
|
||
|
|
// -- premier incrément
|
||
|
|
val_finale = val_initiale + (deltat * bb21) * der_initiale;
|
||
|
|
//// --- debug ---
|
||
|
|
//cout << "\n Algo_edp::Runge_Kutta_step23(.. avant appel "
|
||
|
|
// << "\n val_initiale= ";val_initiale.Affiche();
|
||
|
|
//cout << " der_initiale= "; der_initiale.Affiche();
|
||
|
|
//cout << " val_finale= "; val_finale.Affiche();
|
||
|
|
//cout << "\n t0= "<< t0 << ", deltat= "<< deltat<< ", t0+aa2*deltat= "<< (t0+aa2*deltat);
|
||
|
|
//cout << endl;
|
||
|
|
//// --- fin debug ---
|
||
|
|
// appel de la dérivée pour le premier incrément
|
||
|
|
(instance.*Ptder_fonc) (t0+aa2*deltat,val_finale,k2,erreur); k2 *= deltat;
|
||
|
|
//// --- debug ---
|
||
|
|
//cout << "\n Algo_edp::Runge_Kutta_step23(.. après appel "
|
||
|
|
// << "\n k2= ";k2.Affiche();
|
||
|
|
//cout << endl;
|
||
|
|
//// --- fin debug ---
|
||
|
|
if (erreur) // cas où le calcul n'est pas licite
|
||
|
|
{ estime_erreur.Inita(ConstMath::tresgrand); // on indique une erreur monstrueuse
|
||
|
|
val_finale = val_initiale; // on retourne la valeur initiale
|
||
|
|
return;
|
||
|
|
};
|
||
|
|
// -- deuxième incrément
|
||
|
|
val_finale = val_initiale + (deltat*bb31) * der_initiale + bb32 * k2;
|
||
|
|
// appel de la dérivée pour le deuxième incrément (aa3=1.)
|
||
|
|
(instance.*Ptder_fonc) (t0+deltat,val_finale,k3,erreur); k3 *= deltat;
|
||
|
|
if (erreur) // cas où le calcul n'est pas licite
|
||
|
|
{ estime_erreur.Inita(ConstMath::tresgrand); // on indique une erreur monstrueuse
|
||
|
|
val_finale = val_initiale; // on retourne la valeur initiale
|
||
|
|
return;
|
||
|
|
};
|
||
|
|
|
||
|
|
// somme pondéré des incréments
|
||
|
|
val_finale = val_initiale + (deltat*AB1)*der_initiale+AB2*k2+AB3*k3;
|
||
|
|
// estimation de l'erreur comme différence entre les méthodes du quatrième
|
||
|
|
// et cinquième ordre
|
||
|
|
estime_erreur = (deltat*dA1)*der_initiale+dA2*k2+dA3*k3;
|
||
|
|
// on vérifie que la solution finale est licite
|
||
|
|
(instance.*Ptverif_fonc) (t0+deltat,val_finale,erreur);
|
||
|
|
if (erreur) // cas où le calcul n'est pas licite
|
||
|
|
{ estime_erreur.Inita(ConstMath::tresgrand); // on indique une erreur monstrueuse
|
||
|
|
val_finale = val_initiale; // on retourne la valeur initiale
|
||
|
|
return;
|
||
|
|
};
|
||
|
|
|
||
|
|
//// --- debug ---
|
||
|
|
//cout << "\n Algo_edp::Runge_Kutta_step23(.. à la fin, tout c'est bien passée "
|
||
|
|
// << "\n val_finale= ";val_finale.Affiche();
|
||
|
|
//cout << " estime_erreur= "; estime_erreur.Affiche();
|
||
|
|
//cout << "\n t0+deltat= "<< t0+deltat ;
|
||
|
|
//cout << endl;
|
||
|
|
//// --- fin debug ---
|
||
|
|
}
|
||
|
|
|
||
|
|
// 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
|
||
|
|
template <class T>
|
||
|
|
void Algo_edp::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écup de la dimension
|
||
|
|
int taille = val_initiale.Taille();
|
||
|
|
int erreur=0; // init de l'indicateur d'erreur dans le calcul de la fonction
|
||
|
|
// vérif de la dimmension
|
||
|
|
if (taille != k2.Taille())
|
||
|
|
{ k2.Change_taille(taille);k3.Change_taille(taille);
|
||
|
|
k4.Change_taille(taille),k5.Change_taille(taille);
|
||
|
|
// val_finale.Change_taille(taille);
|
||
|
|
};
|
||
|
|
// -- premier incrément
|
||
|
|
val_finale = val_initiale + (deltat * b_b21) * der_initiale;
|
||
|
|
// appel de la dérivée pour le premier incrément
|
||
|
|
(instance.*Ptder_fonc) (t0+a_a2*deltat,val_finale,k2,erreur); k2 *= deltat;
|
||
|
|
if (erreur) // cas où le calcul n'est pas licite
|
||
|
|
{ estime_erreur.Inita(ConstMath::tresgrand); // on indique une erreur monstrueuse
|
||
|
|
val_finale = val_initiale; // on retourne la valeur initiale
|
||
|
|
return;
|
||
|
|
};
|
||
|
|
// -- deuxième incrément
|
||
|
|
val_finale = val_initiale + (deltat * b_b31) * der_initiale + b_b32* k2;
|
||
|
|
// appel de la dérivée pour le deuxième incrément
|
||
|
|
(instance.*Ptder_fonc) (t0+a_a3*deltat,val_finale,k3,erreur); k3 *= deltat;
|
||
|
|
if (erreur) // cas où le calcul n'est pas licite
|
||
|
|
{ estime_erreur.Inita(ConstMath::tresgrand); // on indique une erreur monstrueuse
|
||
|
|
val_finale = val_initiale; // on retourne la valeur initiale
|
||
|
|
return;
|
||
|
|
};
|
||
|
|
// -- troisième incrément
|
||
|
|
val_finale = val_initiale + (deltat*b_b41)*der_initiale+b_b42*k2+b_b43*k3;
|
||
|
|
// appel de la dérivée pour le troisième incrément
|
||
|
|
(instance.*Ptder_fonc) (t0+a_a4*deltat,val_finale,k4,erreur); k4 *= deltat;
|
||
|
|
if (erreur) // cas où le calcul n'est pas licite
|
||
|
|
{ estime_erreur.Inita(ConstMath::tresgrand); // on indique une erreur monstrueuse
|
||
|
|
val_finale = val_initiale; // on retourne la valeur initiale
|
||
|
|
return;
|
||
|
|
};
|
||
|
|
// -- quatrième incrément
|
||
|
|
val_finale = val_initiale + (deltat*b_b51)*der_initiale+b_b53*k3+b_b54*k4;
|
||
|
|
// appel de la dérivée pour le quatrième incrément
|
||
|
|
(instance.*Ptder_fonc) (t0+a_a5*deltat,val_finale,k5,erreur); k5 *= deltat;
|
||
|
|
if (erreur) // cas où le calcul n'est pas licite
|
||
|
|
{ estime_erreur.Inita(ConstMath::tresgrand); // on indique une erreur monstrueuse
|
||
|
|
val_finale = val_initiale; // on retourne la valeur initiale
|
||
|
|
return;
|
||
|
|
};
|
||
|
|
// somme pondéré des incréments
|
||
|
|
val_finale = val_initiale + (deltat*A_B1)*der_initiale+A_B3*k3+A_B4*k4+A_B5*k5;
|
||
|
|
// on vérifie que la solution finale est licite
|
||
|
|
(instance.*Ptverif_fonc) (t0+deltat,val_finale,erreur);
|
||
|
|
if (erreur) // cas où le calcul n'est pas licite
|
||
|
|
{ estime_erreur.Inita(ConstMath::tresgrand); // on indique une erreur monstrueuse
|
||
|
|
val_finale = val_initiale; // on retourne la valeur initiale
|
||
|
|
return;
|
||
|
|
};
|
||
|
|
// estimation de l'erreur comme différence entre les méthodes du quatrième
|
||
|
|
// et cinquième ordre
|
||
|
|
estime_erreur = (deltat*d_A1)*der_initiale+d_A3*k3+d_A4*k4+d_A5*k5;
|
||
|
|
}
|
||
|
|
|
||
|
|
// 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
|
||
|
|
template <class T>
|
||
|
|
void Algo_edp::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)
|
||
|
|
{// récup de la dimension
|
||
|
|
int taille = val_initiale.Taille();
|
||
|
|
int erreur=0; // init de l'indicateur d'erreur dans le calcul de la fonction
|
||
|
|
// vérif de la dimmension
|
||
|
|
if (taille != k2.Taille())
|
||
|
|
{ k2.Change_taille(taille);k3.Change_taille(taille);
|
||
|
|
k4.Change_taille(taille),k5.Change_taille(taille);
|
||
|
|
k6.Change_taille(taille);//f_inter.Change_taille(taille);
|
||
|
|
};
|
||
|
|
// -- premier incrément
|
||
|
|
val_finale = val_initiale + (deltat * b21) * der_initiale;
|
||
|
|
// appel de la dérivée pour le premier incrément
|
||
|
|
(instance.*Ptder_fonc) (t0+a2*deltat,val_finale,k2,erreur); k2 *= deltat;
|
||
|
|
if (erreur) // cas où le calcul n'est pas licite
|
||
|
|
{ estime_erreur.Inita(ConstMath::tresgrand); // on indique une erreur monstrueuse
|
||
|
|
val_finale = val_initiale; // on retourne la valeur initiale
|
||
|
|
return;
|
||
|
|
};
|
||
|
|
// -- deuxième incrément
|
||
|
|
val_finale = val_initiale + (deltat * b31) * der_initiale + b32* k2;
|
||
|
|
// appel de la dérivée pour le deuxième incrément
|
||
|
|
(instance.*Ptder_fonc) (t0+a3*deltat,val_finale,k3,erreur); k3 *= deltat;
|
||
|
|
if (erreur) // cas où le calcul n'est pas licite
|
||
|
|
{ estime_erreur.Inita(ConstMath::tresgrand); // on indique une erreur monstrueuse
|
||
|
|
val_finale = val_initiale; // on retourne la valeur initiale
|
||
|
|
return;
|
||
|
|
};
|
||
|
|
// -- troisième incrément
|
||
|
|
val_finale = val_initiale + (deltat*b41)*der_initiale+b42*k2+b43*k3;
|
||
|
|
// appel de la dérivée pour le troisième incrément
|
||
|
|
(instance.*Ptder_fonc) (t0+a4*deltat,val_finale,k4,erreur); k4 *= deltat;
|
||
|
|
if (erreur) // cas où le calcul n'est pas licite
|
||
|
|
{ estime_erreur.Inita(ConstMath::tresgrand); // on indique une erreur monstrueuse
|
||
|
|
val_finale = val_initiale; // on retourne la valeur initiale
|
||
|
|
return;
|
||
|
|
};
|
||
|
|
// -- quatrième incrément
|
||
|
|
val_finale = val_initiale + (deltat*b51)*der_initiale+b52*k2+b53*k3+b54*k4;
|
||
|
|
// appel de la dérivée pour le quatrième incrément
|
||
|
|
(instance.*Ptder_fonc) (t0+a5*deltat,val_finale,k5,erreur); k5 *= deltat;
|
||
|
|
if (erreur) // cas où le calcul n'est pas licite
|
||
|
|
{ estime_erreur.Inita(ConstMath::tresgrand); // on indique une erreur monstrueuse
|
||
|
|
val_finale = val_initiale; // on retourne la valeur initiale
|
||
|
|
return;
|
||
|
|
};
|
||
|
|
// -- cinquième incrément
|
||
|
|
val_finale = val_initiale + (deltat*b61)*der_initiale+b62*k2+b63*k3+b64*k4+b65*k5;
|
||
|
|
// appel de la dérivée pour le cinquième incrément
|
||
|
|
(instance.*Ptder_fonc) (t0+a6*deltat,val_finale,k6,erreur); k6 *= deltat;
|
||
|
|
if (erreur) // cas où le calcul n'est pas licite
|
||
|
|
{ estime_erreur.Inita(ConstMath::tresgrand); // on indique une erreur monstrueuse
|
||
|
|
val_finale = val_initiale; // on retourne la valeur initiale
|
||
|
|
return;
|
||
|
|
};
|
||
|
|
// somme pondéré des incréments
|
||
|
|
val_finale = val_initiale + (deltat*c1)*der_initiale+c3*k3+c4*k4+c6*k6;
|
||
|
|
// on vérifie que la solution finale est licite
|
||
|
|
(instance.*Ptverif_fonc) (t0+deltat,val_finale,erreur);
|
||
|
|
if (erreur) // cas où le calcul n'est pas licite
|
||
|
|
{ estime_erreur.Inita(ConstMath::tresgrand); // on indique une erreur monstrueuse
|
||
|
|
val_finale = val_initiale; // on retourne la valeur initiale
|
||
|
|
return;
|
||
|
|
};
|
||
|
|
// estimation de l'erreur comme différence entre les méthodes du quatrième
|
||
|
|
// et cinquième ordre
|
||
|
|
estime_erreur = (deltat*ce1)*der_initiale+ce3*k3+ce4*k4+ce5*k5+ce6*k6;
|
||
|
|
}
|
||
|
|
|
||
|
|
// 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
|
||
|
|
// 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
|
||
|
|
// =9 : tout c'est bien passé, mais la dernière valeur à t, n'a pu être calculée, donc c'est la valeur à t - quelque chose qui est renvoyée
|
||
|
|
// = 0 : erreur inconnue
|
||
|
|
|
||
|
|
template <class T>
|
||
|
|
int Algo_edp::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)
|
||
|
|
{ // -- premières initialisation
|
||
|
|
nb_step=0; // compteur d'incrément
|
||
|
|
// double temps = tfin ; //******* modif 4 septembre:
|
||
|
|
double temps = tdebut;
|
||
|
|
nombreAppelF=0; // compteur global d'appel de Kutta
|
||
|
|
double deltat=tfin-tdebut;
|
||
|
|
erreur_maxi_global = 0.;
|
||
|
|
bool increment_valide=false; // pour savoir si l'exé d'un increment est valide
|
||
|
|
// -- on commence par faire quelques tests pour vérifier le contexte
|
||
|
|
// 1) vérification du niveau de la précision relative demandée
|
||
|
|
static double relerr_min = 2. * ConstMath::eps_machine * remin;
|
||
|
|
|
||
|
|
//// --- debug ---
|
||
|
|
//cout << "\n Algo_edp::Pilotage_kutta( "
|
||
|
|
// << "remin = " << remin << " eps_machine = " << ConstMath::eps_machine << " relerr_min= "
|
||
|
|
// << relerr_min << endl;
|
||
|
|
//
|
||
|
|
//// --- fin debug ---
|
||
|
|
|
||
|
|
|
||
|
|
if (erreurRelative < relerr_min)
|
||
|
|
{ if (permet_affichage > 4 )
|
||
|
|
cout << "\n Erreur: au niveau du pilotage de Runge Kutta, la precision relative demande= "
|
||
|
|
<< erreurRelative << " est inferieur a celle permise= " << relerr_min
|
||
|
|
<< "\n Algo_edp::Pilotage_kutta(..";
|
||
|
|
erreurRelative = relerr_min;
|
||
|
|
return 3; // indication au programme appelant
|
||
|
|
};
|
||
|
|
// -- récup de la dimension
|
||
|
|
int taille = val_initiale.Taille();
|
||
|
|
// vérif de la dimmension
|
||
|
|
if (taille != val_inter.Taille())
|
||
|
|
{ val_inter.Change_taille(taille);val_finale.Change_taille(taille);
|
||
|
|
der_inter.Change_taille(taille);der_finale.Change_taille(taille);
|
||
|
|
};
|
||
|
|
// -- initialisation
|
||
|
|
// -- initialisation avec sauvegarde des valeurs initiales
|
||
|
|
der_inter=der_initiale;
|
||
|
|
val_inter=val_initiale;
|
||
|
|
// def du minimum de step possible (cf. rkf45)
|
||
|
|
double mini_delta_t = 26. * ConstMath::eps_machine * MaX(Dabs(tdebut),Dabs(tfin));
|
||
|
|
// dans cette boucle l'idée est de piloter le pas de temps
|
||
|
|
while (Dabs(temps) < Dabs(tfin))
|
||
|
|
{ increment_valide = true; // par défaut on est optimiste
|
||
|
|
// on appel la routine de calcul de base,
|
||
|
|
// (val_inter, der_inter, temps et deltat : ne sont pas modifiés par le calcul)
|
||
|
|
switch (cas_kutta)
|
||
|
|
{ case 3: // cas kutta23
|
||
|
|
{this->Runge_Kutta_step23(instance,Ptder_fonc,Ptverif_fonc,val_inter,der_inter,temps,deltat
|
||
|
|
,val_finale,estime_erreur);
|
||
|
|
break;
|
||
|
|
};
|
||
|
|
case 4: // cas kutta34
|
||
|
|
{this->Runge_Kutta_step34(instance,Ptder_fonc,Ptverif_fonc,val_inter,der_inter,temps,deltat
|
||
|
|
,val_finale,estime_erreur);
|
||
|
|
break;
|
||
|
|
};
|
||
|
|
case 5: // cas kutta45
|
||
|
|
{this->Runge_Kutta_step45(instance,Ptder_fonc,Ptverif_fonc,val_inter,der_inter,temps,deltat
|
||
|
|
,val_finale,estime_erreur);
|
||
|
|
break;
|
||
|
|
};
|
||
|
|
default:
|
||
|
|
{ if (permet_affichage >4)
|
||
|
|
cout << "\n Erreur: au niveau du pilotage de Runge Kutta, le cas kutta demande= "
|
||
|
|
<< cas_kutta << " est errone ! par de routine de base de ce type (cf. doc) "
|
||
|
|
<< "\n Algo_edp::Pilotage_kutta(..";
|
||
|
|
return 8; // indication au programme appelant
|
||
|
|
};
|
||
|
|
};
|
||
|
|
nombreAppelF += cas_kutta; //5; // on incrémente le compteur d'appels
|
||
|
|
|
||
|
|
// calcul de la nouvelle dérivée initiale pour l'itération éventuelle suivante
|
||
|
|
// on le fait ici, pour vérifier à ce stade que la valeur finale et la dérivée est licite
|
||
|
|
// si ce n'est pas le cas cela veut dire que le calcul précédent n'est finalement pas bon
|
||
|
|
// cela peut-être le cas par exemple lorsque l'évolution est instable (le pas de temps deltat est trop grand)
|
||
|
|
// .. dans le cas où il y a eu une erreur dans le RKij précédent, val_finale est remis à la valeur initiale
|
||
|
|
// c-a-d à val_initiale (représenté par val_inter, dans l'appel de RKij)
|
||
|
|
int erreur_fonc=0; // initialisation
|
||
|
|
// --- debug contrôlé---
|
||
|
|
if (permet_affichage > 6)
|
||
|
|
{cout << "\n Algo_edp::Pilotage_kutta(.. juste après int erreur_fonc=0 "
|
||
|
|
<< "\n le vecteur= ";val_finale.Affiche();
|
||
|
|
cout << "\n l'erreur estime= ";estime_erreur.Affiche();
|
||
|
|
cout << endl;
|
||
|
|
};
|
||
|
|
// --- fin debug ---
|
||
|
|
// calcul de la dérivée finale, car on ne dispose pour l'instant que la la valeur finale
|
||
|
|
(instance.*Ptder_fonc) (temps,val_finale,der_finale,erreur_fonc);
|
||
|
|
// --- debug contrôlé ---
|
||
|
|
if (permet_affichage > 6)
|
||
|
|
{cout << "\n Algo_edp::Pilotage_kutta(.. appel de la fct "
|
||
|
|
<< "\n le vecteur= ";val_finale.Affiche();
|
||
|
|
cout << " sa dérivée= "; der_finale.Affiche();
|
||
|
|
cout << "\n temps= "<< temps <<" et temps finale visee= "<<(temps+deltat);
|
||
|
|
cout << endl;
|
||
|
|
};
|
||
|
|
// --- fin debug ---
|
||
|
|
nombreAppelF++; // incrémente le compteur d'appel de fonction
|
||
|
|
dernierTemps = temps; dernierdeltat=deltat; // init et en même temps sert de sauvegarde
|
||
|
|
|
||
|
|
// ... on vérifie que le nombre d'appel n'exède pas le maxi fixé
|
||
|
|
if (nombreAppelF > nbMaxiAppel)
|
||
|
|
{ if (permet_affichage >4)
|
||
|
|
cout << "\n Erreur au niveau du pilotage de Runge Kutta, le nombre maxi d'appel= "
|
||
|
|
<< nbMaxiAppel << " de la fonction est depasse "
|
||
|
|
<< "\n Algo_edp::Pilotage_kutta(..";
|
||
|
|
if (erreur_fonc) // si erreur dans le dernier calcul, cela signifie que la valeur calculée n'est pas licite
|
||
|
|
// dans ce cas on renvoie la valeur précédente
|
||
|
|
{ dernierTemps -= dernierdeltat;val_finale = val_inter;der_finale = der_inter;}
|
||
|
|
// sinon on renvoie la dernière valeur calculée, qui est licite
|
||
|
|
return 4; // indication au programme appelant
|
||
|
|
};
|
||
|
|
// ... si le dernier calcul a posé pb, la seule solution est de diminuer le pas de temps
|
||
|
|
// dans dernier calcul: on entend le calcul de la dérivée initiale et/ou la dernière résolution
|
||
|
|
// de manière arbitraire et de recommencer le calcul
|
||
|
|
if ((erreur_fonc)|| (estime_erreur.Max_val_abs() > ConstMath::grand))
|
||
|
|
{ deltat /= 2.; // on divise par deux de manière arbitraire
|
||
|
|
//// --- debug ---
|
||
|
|
//cout << "\n Algo_edp::Pilotage_kutta( debut: erreur_fonc= " << erreur_fonc << " estime_erreur.Max_val_abs()= "
|
||
|
|
// << estime_erreur.Max_val_abs() << " deltat= " << deltat;
|
||
|
|
//// --- fin debug ---
|
||
|
|
// on test néanmoins pour savoir si l'on n'est pas en dessous du delta limite inférieur
|
||
|
|
if (Dabs(deltat) < mini_delta_t )
|
||
|
|
{ if (permet_affichage > 4 )
|
||
|
|
cout << "\n Erreur1 au niveau du pilotage de Runge Kutta, le nouvelle increment qu'il faudrait utiliser = "
|
||
|
|
<< deltat << " est plus petit que le minimum autorise: " << mini_delta_t
|
||
|
|
<< "\n Algo_edp::Pilotage_kutta(..";
|
||
|
|
dernierTemps = temps - 2. * deltat;val_finale = val_inter;der_finale = der_inter;
|
||
|
|
return 6; // indication au programme appelant
|
||
|
|
};
|
||
|
|
}
|
||
|
|
else // sinon on peu continuer
|
||
|
|
{
|
||
|
|
// ... on passe en revue toutes les erreurs individuelles de chaque fonction
|
||
|
|
// pour avoir le maximum, comparé à l'erreur admise -> facteur
|
||
|
|
double facteur=0.; // facteur fonction des erreurs: pour contrôler le pas
|
||
|
|
double max_erreur_sur_step_estimee=0;
|
||
|
|
//// --- debug ---
|
||
|
|
//cout << "\n Algo_edp::Pilotage_kutta( ";
|
||
|
|
//// --- fin debug ---
|
||
|
|
|
||
|
|
for (int i=1; i<= taille; i++)
|
||
|
|
{ double erreur_admise = erreurAbsolue+Dabs(val_finale(i))*erreurRelative;
|
||
|
|
// comme on va diviser par erreur_admise, s'il est nulle on l'augmente
|
||
|
|
erreur_admise = MaX(erreur_admise,ConstMath::pasmalpetit);
|
||
|
|
facteur = MaX(facteur,Dabs(estime_erreur(i))/erreur_admise);
|
||
|
|
max_erreur_sur_step_estimee = MaX(max_erreur_sur_step_estimee,Dabs(estime_erreur(i)));
|
||
|
|
//// --- debug ---
|
||
|
|
//cout << "\n estime_erreur(i)= " << estime_erreur(i) << " erreur_admise= "<< erreur_admise
|
||
|
|
// << " max_erreur_sur_step_estimee= "<< max_erreur_sur_step_estimee;
|
||
|
|
//// --- fin debug ---
|
||
|
|
};
|
||
|
|
// on ajuste maintenant le pas en fonction du facteur qui dépend de l'erreur calculée
|
||
|
|
double AbsDeltat=Abs(deltat);
|
||
|
|
if (facteur > 1.)
|
||
|
|
{// signifie que l'erreur du calcul est supérieure à celle admise
|
||
|
|
// on doit donc diminuer le pas et recommencer le calcul
|
||
|
|
double deltat_essai=coef*AbsDeltat*pow(facteur,-0.25); // le delta essai
|
||
|
|
// on limite l'affaiblissement à 1/10 cf rkf45
|
||
|
|
deltat = DSigne(deltat)*MaX(0.1*AbsDeltat,deltat_essai);
|
||
|
|
increment_valide = false;
|
||
|
|
|
||
|
|
//// --- debug ---
|
||
|
|
//cout << "\n Algo_edp::Pilotage_kutta( "
|
||
|
|
// << "facteur = " << facteur << " max_erreur_sur_step_estimee = " << max_erreur_sur_step_estimee << " deltat_essai= "
|
||
|
|
// << deltat_essai << " deltat= "<< deltat << endl;
|
||
|
|
//
|
||
|
|
//// --- fin debug ---
|
||
|
|
|
||
|
|
// on test pour savoir si l'on n'est pas en dessous du delta limite inférieur
|
||
|
|
if (Dabs(deltat) < mini_delta_t )
|
||
|
|
{ if (permet_affichage > 4 )
|
||
|
|
cout << "\n Erreur2 au niveau du pilotage de Runge Kutta, le nouvelle increment qu'il faudrait utiliser = "
|
||
|
|
<< deltat << " est plus petit que le minimum autorise: " << mini_delta_t
|
||
|
|
<< "\n Algo_edp::Pilotage_kutta(..";
|
||
|
|
if (permet_affichage > 5 )
|
||
|
|
cout << "\n Algo_edp::Pilotage_kutta( "
|
||
|
|
<< "facteur = " << facteur << " max_erreur_sur_step_estimee = " << max_erreur_sur_step_estimee << " deltat_essai= "
|
||
|
|
<< deltat_essai << " deltat= "<< deltat << endl;
|
||
|
|
dernierTemps = temps;
|
||
|
|
dernierdeltat=deltat;
|
||
|
|
return 6; // indication au programme appelant
|
||
|
|
};
|
||
|
|
}
|
||
|
|
else // sinon tout est ok: pas de temps, erreur, nombres d'appels et valeurs finales
|
||
|
|
// donc cas ou le calcul à bien fonctionné
|
||
|
|
{// on met à jour le temps
|
||
|
|
temps += deltat; nb_step++;
|
||
|
|
erreur_maxi_global += max_erreur_sur_step_estimee;
|
||
|
|
//// --- debug ---
|
||
|
|
//cout << "\n Algo_edp::Pilotage_kutta( calcul validé: "
|
||
|
|
// << " erreur_maxi_global= "<<erreur_maxi_global << endl;
|
||
|
|
//
|
||
|
|
//// --- fin debug ---
|
||
|
|
if (temps == tfin)
|
||
|
|
{ // dernier calcul pour une fin qui est normale
|
||
|
|
//----------------------------- premier bloc de fin normal ---------------
|
||
|
|
// sinon ------ fin normal du calcul -------
|
||
|
|
// on renvoie la dernière valeur calculée, qui est licite
|
||
|
|
dernierTemps = temps;
|
||
|
|
nombreAppelF++;
|
||
|
|
return 2;
|
||
|
|
//----------------------------- fin du premier bloc de fin normal ---------------
|
||
|
|
};
|
||
|
|
|
||
|
|
// -- si on est trop près de la fin, on extrapole et retour
|
||
|
|
// on utilise une borne défini dans rkf45
|
||
|
|
if (Dabs(tfin-temps) <= 26.* ConstMath::eps_machine * Dabs(deltat))
|
||
|
|
{dernierTemps = temps; // dernier temps licite
|
||
|
|
temps=tfin;
|
||
|
|
der_finale = der_inter; // sauvegarde de la dernière dérivée licite, on se sert de der_finale
|
||
|
|
int erreur1=0;
|
||
|
|
(instance.*Ptder_fonc) (tfin,val_finale,der_inter,erreur1);
|
||
|
|
val_finale += (tfin-temps) * der_inter; // euler
|
||
|
|
if (erreur1)// il y a eu erreur, donc on récupère la dernière dérivée licite
|
||
|
|
der_inter = der_finale;
|
||
|
|
int erreur2=0;
|
||
|
|
(instance.*Ptder_fonc) (tfin,val_finale,der_finale,erreur2);
|
||
|
|
// on regarde s'il y a eu une erreur
|
||
|
|
if (erreur1 || erreur2) // si erreur, cela signifie que la dernière valeur calculée n'est pas licite
|
||
|
|
// dans ce cas on renvoie la valeur précédente
|
||
|
|
{ val_finale = val_inter;der_finale = der_inter;
|
||
|
|
return 9;
|
||
|
|
}
|
||
|
|
//----------------------------- deuxième bloc de fin normal ---------------
|
||
|
|
else
|
||
|
|
{ nombreAppelF +=2; dernierTemps = temps; dernierdeltat=deltat;
|
||
|
|
return 2;
|
||
|
|
};
|
||
|
|
//----------------------------- fin du deuxième bloc de fin normal ---------------
|
||
|
|
};
|
||
|
|
// sinon on peut au contraire augmenter le pas
|
||
|
|
double deltat_essai=coef*AbsDeltat*pow(facteur,-0.2);
|
||
|
|
// on limite l'augmentation à 5 fois le pas initial
|
||
|
|
deltat = DSigne(deltat)*MiN(MiN(5.*AbsDeltat,deltat_essai),Dabs(tfin-temps));
|
||
|
|
//// --- debug ---
|
||
|
|
//cout << "\n Algo_edp::Pilotage_kutta( *** augmentation du pas de temps éventuelle \n "
|
||
|
|
// << " deltat_essai= "
|
||
|
|
// << deltat_essai << " deltat= "<< deltat << " erreur_maxi_global= "<<erreur_maxi_global << endl;
|
||
|
|
//
|
||
|
|
//// --- fin debug ---
|
||
|
|
};
|
||
|
|
};
|
||
|
|
// on initialise l'appel du runge qui va suivre
|
||
|
|
if (increment_valide) // uniquement si l'incrément a été validé
|
||
|
|
{ val_inter = val_finale;
|
||
|
|
der_inter = der_finale;
|
||
|
|
};
|
||
|
|
}; // -- fin du while
|
||
|
|
// pour faire taire le compilo car on ne doit jamais passer ici
|
||
|
|
return 0; // erreur inconnue
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
#endif
|