/*************************************************************************************************************************************/
/* */
/* M E T H O D E S D ' I N T E G R A T I O N N U M E R I Q U E P O U R D E S */
/* S Y S T E M E S D ' E Q U A T I O N S D I F F E R E N T I E L L E S D U P R E M I E R O R D R E : */
/* */
/* */
/* Definition : */
/* */
/* Soit le systeme d'equations differentielles */
/* suivant : */
/* */
/* dx */
/* ---- = F (x,y,z,t) */
/* dt x */
/* */
/* dy */
/* ---- = F (x,y,z,t) */
/* dt y */
/* */
/* dz */
/* ---- = F (x,y,z,t) */
/* dt z */
/* */
/* */
/* Resolution en posant h=dt (voir 'METHODE DE CALCUL NUMERIQUE' de JP. Nougier, page 186) : */
/* */
/* */
/* 1-methode d'Euler (schema explicite du premier ordre) : */
/* */
/* dx = F (x,y,z,t).dt */
/* x */
/* */
/* dy = F (x,y,z,t).dt */
/* y */
/* */
/* dz = F (x,y,z,t).dt */
/* z */
/* */
/* avec 'dt' petit, puis : */
/* */
/* x = x + dx */
/* y = y + dy */
/* z = z + dz */
/* */
/* soit : */
/* */
/* x = x + F (x ,y ,z ,t).h */
/* n+1 n x n n n */
/* */
/* y = y + F (x ,y ,z ,t).h */
/* n+1 n y n n n */
/* */
/* z = z + F (x ,y ,z ,t).h */
/* n+1 n z n n n */
/* */
/* */
/* 2-methode de Runge-Kutta (schema explicite du deuxieme ordre) : */
/* */
/* ^ h */
/* x = x + F (x ,y ,z ,t).- */
/* 1 n x n n n 2 */
/* n+- */
/* 2 */
/* */
/* ^ h */
/* y = y + F (x ,y ,z ,t).- */
/* 1 n y n n n 2 */
/* n+- */
/* 2 */
/* */
/* ^ h */
/* z = z + F (x ,y ,z ,t).- */
/* 1 n z n n n 2 */
/* n+- */
/* 2 */
/* */
/* puis : */
/* */
/* ^ ^ ^ h */
/* x = x + F (x ,y ,z ,t+-).h */
/* n+1 n x 1 1 1 2 */
/* n+- n+- n+- */
/* 2 2 2 */
/* */
/* ^ ^ ^ h */
/* y = y + F (x ,y ,z ,t+-).h */
/* n+1 n y 1 1 1 2 */
/* n+- n+- n+- */
/* 2 2 2 */
/* */
/* ^ ^ ^ h */
/* z = z + F (x ,y ,z ,t+-).h */
/* n+1 n z 1 1 1 2 */
/* n+- n+- n+- */
/* 2 2 2 */
/* */
/* */
/* 3-methode de Runge-Kutta (schema explicite du quatrieme ordre) : */
/* */
/* ^ h */
/* x = x + F (x ,y ,z ,t).- */
/* 1 n x n n n 2 */
/* n+- */
/* 2 */
/* */
/* ^ h */
/* y = y + F (x ,y ,z ,t).- */
/* 1 n y n n n 2 */
/* n+- */
/* 2 */
/* */
/* ^ h */
/* z = z + F (x ,y ,z ,t).- */
/* 1 n z n n n 2 */
/* n+- */
/* 2 */
/* */
/* puis : */
/* */
/* ^ */
/* ^ ^ ^ ^ h h */
/* x = x + F (x ,y ,z ,t+-).- */
/* 1 n x 1 1 1 2 2 */
/* n+- n+- n+- n+- */
/* 2 2 2 2 */
/* */
/* ^ */
/* ^ ^ ^ ^ h h */
/* y = y + F (x ,y ,z ,t+-).- */
/* 1 n y 1 1 1 2 2 */
/* n+- n+- n+- n+- */
/* 2 2 2 2 */
/* */
/* ^ */
/* ^ ^ ^ ^ h h */
/* z = z + F (x ,y ,z ,t+-).- */
/* 1 n z 1 1 1 2 2 */
/* n+- n+- n+- n+- */
/* 2 2 2 2 */
/* */
/* puis : */
/* */
/* ^ ^ ^ */
/* ^ ^ ^ ^ h */
/* x = x + F (x ,y ,z ,t+-).h */
/* n+1 n x 1 1 1 2 */
/* n+- n+- n+- */
/* 2 2 2 */
/* */
/* ^ ^ ^ */
/* ^ ^ ^ ^ h */
/* y = y + F (x ,y ,z ,t+-).h */
/* n+1 n y 1 1 1 2 */
/* n+- n+- n+- */
/* 2 2 2 */
/* */
/* ^ ^ ^ */
/* ^ ^ ^ ^ h */
/* z = z + F (x ,y ,z ,t+-).h */
/* n+1 n z 1 1 1 2 */
/* n+- n+- n+- */
/* 2 2 2 */
/* */
/* puis : */
/* */
/* ^ ^ ^ */
/* ^ ^ ^ h ^ ^ ^ h ^ ^ ^ h */
/* x = x + [F (x ,y ,z ,t) + 2.F (x ,y ,z ,t+-) + 2.F (x ,y ,z ,t+-) + F (x ,y ,z ,t+h)].- */
/* n+1 n x n n n x 1 1 1 2 x 1 1 1 2 x n+1 n+1 n+1 6 */
/* n+- n+- n+- n+- n+- n+- */
/* 2 2 2 2 2 2 */
/* */
/* ^ ^ ^ */
/* ^ ^ ^ h ^ ^ ^ h ^ ^ ^ h */
/* y = y + [F (x ,y ,z ,t) + 2.F (x ,y ,z ,t+-) + 2.F (x ,y ,z ,t+-) + F (x ,y ,z ,t+h)].- */
/* n+1 n y n n n y 1 1 1 2 y 1 1 1 2 y n+1 n+1 n+1 6 */
/* n+- n+- n+- n+- n+- n+- */
/* 2 2 2 2 2 2 */
/* */
/* ^ ^ ^ */
/* ^ ^ ^ h ^ ^ ^ h ^ ^ ^ h */
/* z = z + [F (x ,y ,z ,t) + 2.F (x ,y ,z ,t+-) + 2.F (x ,y ,z ,t+-) + F (x ,y ,z ,t+h)].- */
/* n+1 n z n n n z 1 1 1 2 z 1 1 1 2 z n+1 n+1 n+1 6 */
/* n+- n+- n+- n+- n+- n+- */
/* 2 2 2 2 2 2 */
/* */
/* */
/* Author of '$xrk/integr.1B$vv$I' : */
/* */
/* Jean-Francois Colonna (LACTAMME, 1992??????????). */
/* */
/*************************************************************************************************************************************/
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* D O N N E E S U T I L E S : */
/* */
/*************************************************************************************************************************************/
DEFV(Local,DEFV(Float,INIT(cx,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(Icx,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(cx_t_1,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(cx1,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(cxc1,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(cxc2,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(dcx,FLOT__UNDEF)));
/* Definition de 'x' et 'dx' et des approximations "chapeau" et "chapeau-chapeau" : */
/* */
/* cx_t_1 : x */
/* n */
/* */
/* Icx : valeur intermediaire due a 'INTEGRE(...)' */
/* */
/* ^ */
/* cx1 : x */
/* n+1 */
/* */
/* ^ */
/* cxc1 : x */
/* 1 */
/* n+- */
/* 2 */
/* */
/* ^ */
/* ^ */
/* cxc2 : x */
/* 1 */
/* n+- */
/* 2 */
/* */
DEFV(Local,DEFV(Float,INIT(cy,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(Icy,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(cy_t_1,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(cy1,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(cyc1,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(cyc2,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(dcy,FLOT__UNDEF)));
/* Definition de 'y' et 'dy' et des approximations "chapeau" et "chapeau-chapeau" : */
/* */
/* cy_t_1 : y */
/* n */
/* */
/* Icy : valeur intermediaire due a 'INTEGRE(...)' */
/* */
/* ^ */
/* cy1 : y */
/* n+1 */
/* */
/* ^ */
/* cyc1 : y */
/* 1 */
/* n+- */
/* 2 */
/* */
/* ^ */
/* ^ */
/* cyc2 : y */
/* 1 */
/* n+- */
/* 2 */
/* */
DEFV(Local,DEFV(Float,INIT(cz,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(Icz,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(cz_t_1,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(cz1,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(czc1,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(czc2,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(dcz,FLOT__UNDEF)));
/* Definition de 'z' et 'dz' et des approximations "chapeau" et "chapeau-chapeau" : */
/* */
/* cy_t_1 : y */
/* n */
/* */
/* Icz : valeur intermediaire due a 'INTEGRE(...)' */
/* */
/* ^ */
/* cz1 : z */
/* n+1 */
/* */
/* ^ */
/* czc1 : z */
/* 1 */
/* n+- */
/* 2 */
/* */
/* ^ */
/* ^ */
/* czc2 : z */
/* 1 */
/* n+- */
/* 2 */
/* */
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* F O R M U L E S D E B A S E : */
/* */
/*************************************************************************************************************************************/
#define INTEGRE(v2,v1,dv,fonction,h) \
Bblock \
EGAL(dv,MUL2(fonction,h)); \
/* Calcul de la differentielle : */ \
/* */ \
/* dv = fonction.h */ \
/* */ \
\
Test(aIFID(v1,v2)) \
Bblock \
PRINT_ERREUR("une coordonnee va etre mise a jour alors que son ancienne valeur peut etre encore utile"); \
/* En effet, lors de l'integration qui va suivre, on calcule : */ \
/* */ \
/* v2 = v1 + dv */ \
/* */ \
/* or, lorsque 'v1' et 'v2' sont identiques, il y a un risque tres fort que l'ancienne */ \
/* valeur 'v1' soit encore utile par la suite (dans l'evaluation des fonctions F(x,y,z,t)), */ \
/* il faut donc prevoir une valeur intermediaire... */ \
Eblock \
ATes \
Bblock \
Eblock \
ETes \
\
EGAL(v2,ADD2(v1,dv)); \
/* Puis integration : */ \
/* */ \
/* v2 = v1 + dv */ \
/* */ \
Eblock \
/* Procedure d'integration... */
#define MISE_A_JOUR(v12,v11,v22,v21,v32,v31) \
Bblock \
EGAL(v12,v11); \
EGAL(v22,v21); \
EGAL(v32,v31); \
Eblock \
/* Procedure de mise a jour a la suite de 'INTEGRE(...)' lorsque l'on est dans le cas ou */ \
/* les variables 'v1' et 'v2' sont identiques... */
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* C H O I X D E L ' O R D R E D E L A M E T H O D E D ' I N T E G R A T I O N : */
/* */
/*************************************************************************************************************************************/
#define METHODE_D_EULER_D_ORDRE_1 \
UN
#define METHODE_DE_RUNGE_KUTTA_D_ORDRE_2 \
DEUX
#define METHODE_DE_RUNGE_KUTTA_D_ORDRE_4 \
QUATRE
#TestADef ORDRE_DE_LA_METHODE_D_INTEGRATION \
METHODE_D_EULER_D_ORDRE_1
DEFV(Local,DEFV(Int,INIT(ordre_de_la_methode_d_integration,ORDRE_DE_LA_METHODE_D_INTEGRATION)));
/* Choix de la methode d'integration parmi les trois proposees ; les valeurs possibles sont */
/* '1', '2' et '4' (elles definissent l'ordre du schema). Cet indicateur a ete introduit */
/* ici le 20070814110146, alors qu'anterieurement a cette date il figurait explicitement */
/* dans les programmes suivants : */
/* */
/* $xrk/SolConnue.11$K */
/* $xrk/SolConnue.21$K */
/* $xrk/fluide_2D.11$K */
/* $xrk/lorenz.11$K */
/* $xrk/lorenz.21$K */
/* $xrk/rdn_walk.52$K */
/* */
/* mais il est plus logique de le mettre ici... */
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* F O R M U L E S G E N E R A L E S D ' I N T E G R A T I O N : */
/* */
/*************************************************************************************************************************************/
#define gINTEGRATION_D_UN_SYSTEME_D_EQUATIONS_DIFFERENTIELLES_O1(cx,cy,cz,t,h) \
Bblock \
MISE_A_JOUR(cx_t_1,cx,cy_t_1,cy,cz_t_1,cz); \
/* Sauvegarde de l'etat anterieur {cx_t_1,cy_t_1,cz_t_1}. */ \
\
Choi(ordre_de_la_methode_d_integration) \
Bblock \
Ca1e(METHODE_D_EULER_D_ORDRE_1) \
Bblock \
INTEGRE(Icx,cx,dcx,Fx(cx,cy,cz,t),FRA1(h)); \
INTEGRE(Icy,cy,dcy,Fy(cx,cy,cz,t),FRA1(h)); \
INTEGRE(Icz,cz,dcz,Fz(cx,cy,cz,t),FRA1(h)); \
/* Integration par la methode d'Euler (ordre 1) : */ \
/* */ \
/* x = x + F (x ,y ,z ,t).h */ \
/* n+1 n x n n n */ \
/* */ \
/* y = y + F (x ,y ,z ,t).h */ \
/* n+1 n y n n n */ \
/* */ \
/* z = z + F (x ,y ,z ,t).h */ \
/* n+1 n z n n n */ \
/* */ \
Eblock \
ECa1 \
\
Ca1e(METHODE_DE_RUNGE_KUTTA_D_ORDRE_2) \
Bblock \
INTEGRE(cxc1,cx,dcx,Fx(cx,cy,cz,t),FRA2(h)); \
INTEGRE(cyc1,cy,dcy,Fy(cx,cy,cz,t),FRA2(h)); \
INTEGRE(czc1,cz,dcz,Fz(cx,cy,cz,t),FRA2(h)); \
/* Calcul des approximations "chapeau" centrees : */ \
/* */ \
/* ^ h */ \
/* x = x + F (x ,y ,z ,t).- */ \
/* 1 n x n n n 2 */ \
/* n+- */ \
/* 2 */ \
/* */ \
/* ^ h */ \
/* y = y + F (x ,y ,z ,t).- */ \
/* 1 n y n n n 2 */ \
/* n+- */ \
/* 2 */ \
/* */ \
/* ^ h */ \
/* z = z + F (x ,y ,z ,t).- */ \
/* 1 n z n n n 2 */ \
/* n+- */ \
/* 2 */ \
/* */ \
INTEGRE(Icx,cx,dcx,Fx(cxc1,cyc1,czc1,ADD2(t,FRA2(h))),FRA1(h)); \
INTEGRE(Icy,cy,dcy,Fy(cxc1,cyc1,czc1,ADD2(t,FRA2(h))),FRA1(h)); \
INTEGRE(Icz,cz,dcz,Fz(cxc1,cyc1,czc1,ADD2(t,FRA2(h))),FRA1(h)); \
/* Integration par la methode de Runge-Kutta (ordre 2) : */ \
/* */ \
/* ^ ^ ^ h */ \
/* x = x + F (x ,y ,z ,t+-).h */ \
/* n+1 n x 1 1 1 2 */ \
/* n+- n+- n+- */ \
/* 2 2 2 */ \
/* */ \
/* ^ ^ ^ h */ \
/* y = y + F (x ,y ,z ,t+-).h */ \
/* n+1 n y 1 1 1 2 */ \
/* n+- n+- n+- */ \
/* 2 2 2 */ \
/* */ \
/* ^ ^ ^ h */ \
/* z = z + F (x ,y ,z ,t+-).h */ \
/* n+1 n z 1 1 1 2 */ \
/* n+- n+- n+- */ \
/* 2 2 2 */ \
/* */ \
Eblock \
ECa1 \
\
Ca1e(METHODE_DE_RUNGE_KUTTA_D_ORDRE_4) \
Bblock \
DEFV(Float,INIT(Fx1,Fx(cx,cy,cz,t))); \
DEFV(Float,INIT(Fx2,FLOT__UNDEF)); \
DEFV(Float,INIT(Fx3,FLOT__UNDEF)); \
DEFV(Float,INIT(Fx4,FLOT__UNDEF)); \
/* Pre-calcul de certaines valeurs de la fonction 'Fx(...)'. Cette optimisation a ete */ \
/* rendue necessaire dans 'v $xrk/rdn_walk.52$K F.fonction.coordonnees.' ou les fonctions */ \
/* 'F(...)' peuvent etre tres "lourdes" a calculer lorsqu'il y a beaucoup de corps... */ \
DEFV(Float,INIT(Fy1,Fy(cx,cy,cz,t))); \
DEFV(Float,INIT(Fy2,FLOT__UNDEF)); \
DEFV(Float,INIT(Fy3,FLOT__UNDEF)); \
DEFV(Float,INIT(Fy4,FLOT__UNDEF)); \
/* Pre-calcul de certaines valeurs de la fonction 'Fy(...)'. Cette optimisation a ete */ \
/* rendue necessaire dans 'v $xrk/rdn_walk.52$K F.fonction.coordonnees.' ou les fonctions */ \
/* 'F(...)' peuvent etre tres "lourdes" a calculer lorsqu'il y a beaucoup de corps... */ \
DEFV(Float,INIT(Fz1,Fz(cx,cy,cz,t))); \
DEFV(Float,INIT(Fz2,FLOT__UNDEF)); \
DEFV(Float,INIT(Fz3,FLOT__UNDEF)); \
DEFV(Float,INIT(Fz4,FLOT__UNDEF)); \
/* Pre-calcul de certaines valeurs de la fonction 'Fz(...)'. Cette optimisation a ete */ \
/* rendue necessaire dans 'v $xrk/rdn_walk.52$K F.fonction.coordonnees.' ou les fonctions */ \
/* 'F(...)' peuvent etre tres "lourdes" a calculer lorsqu'il y a beaucoup de corps... */ \
\
INTEGRE(cxc1,cx,dcx,Fx1,FRA2(h)); \
INTEGRE(cyc1,cy,dcy,Fy1,FRA2(h)); \
INTEGRE(czc1,cz,dcz,Fz1,FRA2(h)); \
/* Calcul des approximations "chapeau" centrees. */ \
/* */ \
/* ^ h */ \
/* x = x + F (x ,y ,z ,t).- */ \
/* 1 n x n n n 2 */ \
/* n+- */ \
/* 2 */ \
/* */ \
/* ^ h */ \
/* y = y + F (x ,y ,z ,t).- */ \
/* 1 n y n n n 2 */ \
/* n+- */ \
/* 2 */ \
/* */ \
/* ^ h */ \
/* z = z + F (x ,y ,z ,t).- */ \
/* 1 n z n n n 2 */ \
/* n+- */ \
/* 2 */ \
/* */ \
EGAL(Fx2,Fx(cxc1,cyc1,czc1,ADD2(t,FRA2(h)))); \
EGAL(Fy2,Fy(cxc1,cyc1,czc1,ADD2(t,FRA2(h)))); \
EGAL(Fz2,Fz(cxc1,cyc1,czc1,ADD2(t,FRA2(h)))); \
/* Pre-calcul de certaines valeurs des fonction 'F?(...)'. */ \
\
INTEGRE(cxc2,cx,dcx,Fx2,FRA2(h)); \
INTEGRE(cyc2,cy,dcy,Fy2,FRA2(h)); \
INTEGRE(czc2,cz,dcz,Fz2,FRA2(h)); \
/* Calcul des approximations "chapeau-chapeau" centrees : */ \
/* */ \
/* ^ */ \
/* ^ ^ ^ ^ h h */ \
/* x = x + F (x ,y ,z ,t+-).- */ \
/* 1 n x 1 1 1 2 2 */ \
/* n+- n+- n+- n+- */ \
/* 2 2 2 2 */ \
/* */ \
/* ^ */ \
/* ^ ^ ^ ^ h h */ \
/* y = y + F (x ,y ,z ,t+-).- */ \
/* 1 n y 1 1 1 2 2 */ \
/* n+- n+- n+- n+- */ \
/* 2 2 2 2 */ \
/* */ \
/* ^ */ \
/* ^ ^ ^ ^ h h */ \
/* z = z + F (x ,y ,z ,t+-).- */ \
/* 1 n z 1 1 1 2 2 */ \
/* n+- n+- n+- n+- */ \
/* 2 2 2 2 */ \
/* */ \
EGAL(Fx3,Fx(cxc2,cyc2,czc2,ADD2(t,FRA2(h)))); \
EGAL(Fy3,Fy(cxc2,cyc2,czc2,ADD2(t,FRA2(h)))); \
EGAL(Fz3,Fz(cxc2,cyc2,czc2,ADD2(t,FRA2(h)))); \
/* Pre-calcul de certaines valeurs des fonction 'F?(...)'. */ \
\
INTEGRE(cx1,cx,dcx,Fx3,FRA1(h)); \
INTEGRE(cy1,cy,dcy,Fy3,FRA1(h)); \
INTEGRE(cz1,cz,dcz,Fz3,FRA1(h)); \
/* Calcul des approximations "chapeau" : */ \
/* */ \
/* ^ ^ ^ */ \
/* ^ ^ ^ ^ h */ \
/* x = x + F (x ,y ,z ,t+-).h */ \
/* n+1 n x 1 1 1 2 */ \
/* n+- n+- n+- */ \
/* 2 2 2 */ \
/* */ \
/* ^ ^ ^ */ \
/* ^ ^ ^ ^ h */ \
/* y = y + F (x ,y ,z ,t+-).h */ \
/* n+1 n y 1 1 1 2 */ \
/* n+- n+- n+- */ \
/* 2 2 2 */ \
/* */ \
/* ^ ^ ^ */ \
/* ^ ^ ^ ^ h */ \
/* z = z + F (x ,y ,z ,t+-).h */ \
/* n+1 n z 1 1 1 2 */ \
/* n+- n+- n+- */ \
/* 2 2 2 */ \
/* */ \
EGAL(Fx4,Fx(cx1,cy1,cz1,ADD2(t,h))); \
EGAL(Fy4,Fy(cx1,cy1,cz1,ADD2(t,h))); \
EGAL(Fz4,Fz(cx1,cy1,cz1,ADD2(t,h))); \
/* Pre-calcul de certaines valeurs des fonction 'F?(...)'. */ \
\
INTEGRE(Icx \
,cx \
,dcx \
,ADD4(GRO1(Fx1) \
,GRO2(Fx2) \
,GRO2(Fx3) \
,GRO1(Fx4) \
) \
,FRA6(h) \
); \
INTEGRE(Icy \
,cy \
,dcy \
,ADD4(GRO1(Fy1) \
,GRO2(Fy2) \
,GRO2(Fy3) \
,GRO1(Fy4) \
) \
,FRA6(h) \
); \
INTEGRE(Icz \
,cz \
,dcz \
,ADD4(GRO1(Fz1) \
,GRO2(Fz2) \
,GRO2(Fz3) \
,GRO1(Fz4) \
) \
,FRA6(h) \
); \
/* Integration par la methode de Runge-Kutta (ordre 4) : */ \
/* */ \
/* h */ \
/* x = x + [...].- */ \
/* n+1 n 6 */ \
/* */ \
/* h */ \
/* y = y + [...].- */ \
/* n+1 n 6 */ \
/* */ \
/* h */ \
/* z = z + [...].- */ \
/* n+1 n 6 */ \
/* */ \
Eblock \
ECa1 \
\
Defo \
Bblock \
PRINT_ERREUR("la methode d'integration n'est pas reconnue"); \
Eblock \
EDef \
Eblock \
ECho \
\
MISE_A_JOUR(cx,Icx,cy,Icy,cz,Icz); \
/* Mise a jour du nouvel etat courant {cx,cy,cz}. */ \
Eblock \
/* Integration du systeme d'equations differentielles suivant une certaine methode pour */ \
/* des variables quelconques {cx,cy,cz}. */
#define INTEGRATION_D_UN_SYSTEME_D_EQUATIONS_DIFFERENTIELLES_O1(t,h) \
Bblock \
gINTEGRATION_D_UN_SYSTEME_D_EQUATIONS_DIFFERENTIELLES_O1(cx,cy,cz,t,h); \
Eblock \
/* Integration du systeme d'equations differentielles pour les variables implicites */ \
/* {cx,cy,cz}. */
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* I N T E G R A T I O N A D A P T A T I V E : */
/* */
/*************************************************************************************************************************************/
#define FAIRE_DE_L_INTEGRATION_ADAPTATIVE \
FAUX
DEFV(Local,DEFV(Logical,INIT(faire_de_l_integration_adaptative,FAIRE_DE_L_INTEGRATION_ADAPTATIVE)));
/* Cet indicateur introduit le 20070814113854 permet de faire de l'integration adaptative, */
/* c'est-a-dire une integration ou le pas de temps est adapte en permanence de facon a */
/* faire que nouvel etat courant soit aussi proche de l'etat anterieur que desire... */
/* */
/* ATTENTION : on notera que dans ce cas, le coloriage se faisant bien souvent a l'aide */
/* de {dcx,dcy,dcz} ('v $xrk/lorenz.11$K memorisation_1_point_07'), ces differentielles */
/* {dcx,dcy,dcz} etant calculee par un produit de la fonction par le pas de temps courant */
/* (dans la procedure 'INTEGRE(...)' ci-dessus), celles-ci ne varieront donc que peu d'un */
/* point au suivant et les couleurs seront donc moins variees. Le 20070814175129, je note */
/* qu'evidemment pour ameliorer cela, il suffit d'utiliser l'option : */
/* */
/* extrema_differentielles_veritables=VRAI */
/* */
/* pour retrouver une riche gamme de couleurs... */
#TestADef SEUIL_DE_LA_DISTANCE_LORS_DE_L_INTEGRATION_ADAPTATIVE \
F_INFINI \
/* Le 20070815095148 ce 'define' est devenu un 'TestADef' et sa valeur est passe de 'FU' */ \
/* a 'F_INFINI'... */
DEFV(Local,DEFV(Float,INIT(seuil_de_la_distance_lors_de_l_integration_adaptative
,SEUIL_DE_LA_DISTANCE_LORS_DE_L_INTEGRATION_ADAPTATIVE
)
)
);
/* Seuil destine a declencher l'adaptation du pas de temps (introduit le 20070814113854). */
#define FACTEUR_DE_REDUCTION_DU_PAS_DE_TEMPS \
GRO11(FRA10(FU))
DEFV(Local,DEFV(Float,INIT(facteur_de_reduction_du_pas_de_temps,FACTEUR_DE_REDUCTION_DU_PAS_DE_TEMPS)));
/* Facteur de reduction du pas de temps (introduit le 20070814114758). On notera qu'il est */
/* imperatif que ce facteur soit evidemment plus grand que 1, tout en etant proche de 1 et */
/* ce afin d'eviter des "discontinuites" visuelles entre des zones calculees pour un certain */
/* pas de temps 'dt' et d'autres zones calculees avec ce meme 'dt' divise par le facteur... */
#define INTEGRATION_ADAPTATIVE_D_UN_SYSTEME_D_EQUATIONS_DIFFERENTIELLES_O1(t,h) \
Bblock \
Test(IL_NE_FAUT_PAS(faire_de_l_integration_adaptative)) \
Bblock \
INTEGRATION_D_UN_SYSTEME_D_EQUATIONS_DIFFERENTIELLES_O1(t,h); \
/* Cas de l'integration non adaptative... */ \
Eblock \
ATes \
Bblock \
DEFV(Logical,INIT(adapter_h,VRAI)); \
DEFV(Float,INIT(h_courant,h)); \
/* Donnees de controle de l'adaptation du pas de temps 'h'... */ \
DEFV(Float,INIT(sauvegarde_de_cx,cx)); \
DEFV(Float,INIT(sauvegarde_de_cy,cy)); \
DEFV(Float,INIT(sauvegarde_de_cz,cz)); \
/* Sauvegarde de l'etat anterieur {cx,cy,cz}... */ \
\
\
Tant(IL_FAUT(adapter_h)) \
Bblock \
INTEGRATION_D_UN_SYSTEME_D_EQUATIONS_DIFFERENTIELLES_O1(t,h_courant); \
/* Integration... */ \
\
Test(IFLT(RdisF3D(cx_t_1,cy_t_1,cz_t_1,cx,cy,cz),seuil_de_la_distance_lors_de_l_integration_adaptative)) \
/* ATTENTION : je note le 20070819104010 qu'en toute logique ce n'est pas la distance */ \
/* euclidienne 'RdisF3D(...)' qu'il faudrait utiliser, mais la distance mesuree (par */ \
/* integration) le long de la trajectoire en cours de calcul ; mais malheureusement, */ \
/* cela serait tres complique. Pour des courbes integrales du type de 'v $xiirk/LORE.K1' */ \
/* cela peut etre un probleme car les points P(t-1) et P(t) peuvent etre tres proches dans */ \
/* l'espace R^3 (distance calculee par 'RdisF3D(...)'), alors qu'ils sont en realite tres */ \
/* eloignes sur l'attracteur lui-meme... */ \
Bblock \
EGAL(adapter_h,FAUX); \
/* Lorsque la distance entre {cx_t_1,cy_t_1,cz_t_1} et {cx,cy,cz} est inferieure au seuil, */ \
/* on arrete l'adpatation... */ \
Eblock \
ATes \
Bblock \
EGAL(cx,sauvegarde_de_cx); \
EGAL(cy,sauvegarde_de_cy); \
EGAL(cz,sauvegarde_de_cz); \
/* Lorsque la distance entre {cx_t_1,cy_t_1,cz_t_1} et {cx,cy,cz} est superieure au seuil, */ \
/* on revient en arriere, */ \
EGAL(h_courant,DIVI(h_courant,facteur_de_reduction_du_pas_de_temps)); \
/* Et on reduit le pas de temps... */ \
Eblock \
ETes \
Eblock \
ETan \
Eblock \
ETes \
Eblock \
/* Integration du systeme d'equations differentielles pour les variables implicites */ \
/* {cx,cy,cz} avec adaptation du pas de temps 'h' de facon a rendre le nouvel etat */ \
/* courant {cx,cy,cz} aussi proche que desire de l'etat anterieur. Ceci fut introduit */ \
/* le 20070814104414... */