probleme avec runge kutta dimension 4

Version imprimable

13/11/2006, 23h35
fab13

probleme avec runge kutta dimension 4

bonsoir,

j'utilise la méthode numérique runge kutta de dimension 4 pour résoudre un système de 4 équations différentielles composées chacune de l'addition d'une dérivée seconde d'une des 4 variables de ce système avec des dérivées premieres de ces memes variables. Le but est de modéliser la trajectoire d'un rayon lumineux partant d'une galaxie (géodésiques) dans un univers euclidien (k=0) en expansion.
La métrique de friedmann walker robertson sert à calculer les coeficients (symboles de christoffel) présents dans chacune des équations différentielles.

Ce rayon lumineux part d'une galaxie repérée par r, theta et phi. la 4ème variable est ct (c la célérité de la lumière et t le temps). Chacune de ces 4 variables est paramétré par une variable arbitraire s dans le système d'équations différentielles. les conditions initiales sont déterminantes sur l'interprétation à faire. j'ai choisi : ct=2c/(3Ho) qui vaut à peu près 3000 Mégaparsec, r=dist (distance de la galaxie d'où part le rayon ), theta=pi/2, phi=pi/2 et pour les dérivées: d(ct)/ds=1, d(r)/ds=-1, d(theta)/ds=0 et d(phi)/ds=0.

ceci permet d'avoir d(r)/d(ct)=-1 pour etre sur que le rayon se dirige bien vers notre galaxie (r diminue par rapport au temps). r et ct sont des variables réduites dans le système différentiel, elles s'expriment en Mégaparsec (1 parsec =3.26 année lumière).

Le probleme est le suivant, avec les conditions initiales ci-dessus, le rayon arrive bien dans notre galaxie mais apres, au lieu de s'éloigner et que r se mette à augmenter, r devient négatif et je ne comprends pas pourquoi car en coordonnées polaires, r est tout le temps positif. Par contre, si je change la condition initiale d(theta)/ds (par exemple au lieu de 0, mettre d(theta)/ds=0.0001), r diminue jusqu'à une faible distance de notre galaxie puis augmente, ce qui est le résultat attendu car on a introduit une légère déviation dans le départ du rayon lumineux.
Alors voilà ma question, est-ce un probleme de divergence lié à la méthode runge kutta qui fait que r devient négatif apres avoir dépassé notre galaxie ?

je joins les deux graphiques qui montrent la différence, en abcisse le temps local de notre galaxie exprimé en Megaparsec (ct), en ordonnée la distance de la galaxie qui s'éloigne (trait plein) et celle du rayon lumineux (pointillé).
sur geo_1.pdf, d(theta)/ds=0 et r devient négatif et sur geo_2.pdf, d(theta)/ds=0.000001 et r, apres une diminution, se met à augmenter.
j'ai pris r0=10000 Mégaparsec pour la distance de départ de la galaxie qui nous envoie le rayon lumineux.

Je vous remercie d'avance pour l'aide que vous pourrez m'apporter.
14/11/2006, 20h12
j.p.mignot

Pouriez vous rappeler + en détail les équations utilisées?
A prioris si Theta, Phi dT/dt et dP/dt =0 alors la symétrie du système devait
supprimer toute possibilité de s'éloigner du rayon alors rien d'étonant à ce que T(t)=P(t)=0 & que r(t) varie r->-r ne serait alors qu'une inversion de phase ( passe du photon qui s'approche à la situation où il s'éloigne).
Je serais donc vivement intéressé à voir les équadif utilisées!

équations différentielles utilisées avec runge kutta dimension 4

voici les 4 équations différentielles utilisées, le programme est écrit en langage C:

r représente la dérivée seconde de chacune des 4 variables, dérivée par rapport au paramètre arbitraire s.

a=ct
ap=d(ct)/ds
b=r
bp=dr/ds
c=theta
cp=d(theta)/ds
d=phi
dp=d(phi)/ds

Ho représente la constante de Hubble et cv la vitesse de la lumière.

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
 
long double F(long double a,long double ap,long double b,long double bp,long double c,long double cp,long double d,long double dp)
{ long double r;
  long double R0=1;
  long double Ho=65000; 
  long double cv=3*powl(10,8);
 
/* r=d^2(ct)/ds^2 */
r=-((2./3)*R0*R0*powl((3*Ho)/(2*cv),4./3)*powl((a),1./3))*(powl(bp,2)+powl(b,2)*powl(cp,2)+powl(b,2)*powl(sinl(c),2)*powl(dp,2));
 
 
  return r;
}
 
long double G(long double a,long double ap,long double b,long double bp,long double c,long double cp,long double d,long double dp)
{ long double r;
  long double Ho=65000;
  long double cv=3*powl(10,8);
 
/* r=d^2(r)/ds^2 */
 
  r=-((4./(3*(a)))*ap*bp)+b*(powl(cp,2)+powl(sinl(c),2)*powl(dp,2));
 
  return r;
}
 
long double H(long double a,long double ap,long double b,long double bp,long double c,long double cp,long double d,long double dp)
{ long double r;
  long double Ho=65000;
  long double cv=3*powl(10,8);
 
 
/* r=d^2(theta)/ds^2 */
 
   r=-(2./b)*bp*cp-(4./(3*(a)))*ap*cp+(sinl(2*c)*powl(dp,2))/2;
 
 
  return r;
}
 
long double I(long double a,long double ap,long double b,long double bp,long double c,long double cp,long double d,long double dp)
{ long double r;
  long double Ho=65000;
  long double cv=3*powl(10,8);
 
/* r=d^2(phi)/ds^2 */
 
  r=-2*dp*((2./(3*(a)))*ap+(1./tanl(c))*cp+(1./b)*bp);
  return r;
}

le système est invariant selon phi.

les conditions initiales sont :

ct=(2*c)/Ho
r=10000 Mpc
theta=pi/2
phi=pi/2

d(ct)/ds=1
dr/ds=-1
d(theta)/ds=0
d(phi)/ds=0

avec ces conditions, r diminue et deviens négatif.

le pas choisi est égal à 0.1 dans l'incrémentation de la méthode runge kutta.

équations différentielles utilisées avec runge kutta dimension 4

voici les 4 équations différentielles utilisées, le programme est écrit en langage C:

la variable r (return r) dans ces 4 fonctions représente la dérivée seconde de chacune des 4 variables, dérivée par rapport au paramètre arbitraire s.

a=ct
ap=d(ct)/ds
b=r
bp=dr/ds
c=theta
cp=d(theta)/ds
d=phi
dp=d(phi)/ds

Ho représente la constante de Hubble et cv la vitesse de la lumière.

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
 
long double F(long double a,long double ap,long double b,long double bp,long double c,long double cp,long double d,long double dp)
{ long double r;
  long double R0=1;
  long double Ho=65000; 
  long double cv=3*powl(10,8);
 
/* r=d^2(ct)/ds^2 */
r=-((2./3)*R0*R0*powl((3*Ho)/(2*cv),4./3)*powl((a),1./3))*(powl(bp,2)+powl(b,2)*powl(cp,2)+powl(b,2)*powl(sinl(c),2)*powl(dp,2));
 
 
  return r;
}
 
long double G(long double a,long double ap,long double b,long double bp,long double c,long double cp,long double d,long double dp)
{ long double r;
  long double Ho=65000;
  long double cv=3*powl(10,8);
 
/* r=d^2(r)/ds^2 */
 
  r=-((4./(3*(a)))*ap*bp)+b*(powl(cp,2)+powl(sinl(c),2)*powl(dp,2));
 
  return r;
}
 
long double H(long double a,long double ap,long double b,long double bp,long double c,long double cp,long double d,long double dp)
{ long double r;
  long double Ho=65000;
  long double cv=3*powl(10,8);
 
 
/* r=d^2(theta)/ds^2 */
 
   r=-(2./b)*bp*cp-(4./(3*(a)))*ap*cp+(sinl(2*c)*powl(dp,2))/2;
 
 
  return r;
}
 
long double I(long double a,long double ap,long double b,long double bp,long double c,long double cp,long double d,long double dp)
{ long double r;
  long double Ho=65000;
  long double cv=3*powl(10,8);
 
/* r=d^2(phi)/ds^2 */
 
  r=-2*dp*((2./(3*(a)))*ap+(1./tanl(c))*cp+(1./b)*bp);
  return r;
}