Discussion :

[fab13]

bonsoir, j'ai développé en langage C il y a quelques temps une application visant à déterminer la trajectoire d'un rayon lumineux dans un espace euclidien en expansion (métrique de Friemann Walker Robertson). J'ai utilisé la méthode de runge kutta à l'ordre 4 pour résoudre numériquement un système de 4 équations différentielles du second ordre ( on peut donc se ramener à 8 équations du 1er ordre dans rk4). Les inconnues sont : le temps (c*t en fait), rho, theta et phi (elles sont toutes les 4 fonction d'un paramètre s dans ce système d'ode ). J'obtiens la figure ( cf geo2.pdf ) qui représente la trajectoire d'un rayon lumineux (en abcisse le temps exprimé en Megaparsec et en ordonnée la distance de la galaxie qui s'éloigne (trait plein) et celle du rayon lumineux (pointillé)). le rayon part d'une distance de 3000 Megaparsec.

J'ai essayé de retrouver le même graphique en langage matlab avec la fonction ode45 mais il n'y a que la première partie du graphique qui ressemble à celui obtenu en C (quand le rayon se rapproche de notre galaxie dr/dt <0). la deuxième partie (quand il s'éloigne dr/dt >0) ne semble pas etre correcte ( cf geo3.pdf ). En effet, la courbe de rho qui doit rester positif n'est pas symétrique au prolongement "virtuel" de la première partie de la courbe. Tout ceci ne se produit que si je prends la condition initiale d(theta)/ds différent de 0 . Avec d(theta)/ds=0, je retrouve le graphique n°1. Je signale que dans le cas où je prends la C.I d(theta)/ds différent de 0, j'obtiens malgré tout le déphasage de pi de l'inconnue theta ( cf geo4.pdf ).

Y'a t-il un problème d'intégration dans ce système concernant l'inconnue theta ? Comment se fait-il que la deuxième partie de la courbe soit fausse et comment y remédier ? (utiliser dde23 ? )

Je vous remercie par avance.

ps: le contexte du sujet est expliqué aussi ici :http://www.developpez.net/forums/d23...dimension-4-a/

[FR119492]

Salut!
La première idée qui me vient à l'esprit est que ton système différentiel pourrait être raide (stiff). Dans ce cas, deux méthodes d'intégration peuvent donner des résultats différents. Reste à savoir lequel est faux.
Jean-Marc Blanc

[fab13]

Les résultats que j'obtiens par le programme en langage C sont physiquement corrects. J'y ai implémenté la méthode de runge-kutta à l'ordre 4, ce qui est beaucoup plus long qu'en utilisant le solver ode45 ou ode15 dans matlab qui est censé donner les mêmes solutions. Si mon système est "stiff", quelles solutions me reste t-il pour corriger les erreurs données avec ode45 ?

Voici le code matlab représentant le système d'ode à résoudre numériquement
avec la condition initiale d(theta)/ds différent de 0 :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
 
function geodesiques2
format long;
 
R0=1;
H0=65000;
cv=3*10^(8);
ctk=(2*cv)/(3*H0);
w1=ctk;
x1=3000;
y1=pi/2;
z1=pi/2;
 
xp=-1;
yp=0.001; % d(theta)/ds différent de 0
zp=0;
 
 
wp=sqrt((R0*R0*(3*H0)/(2*cv)*w1)^(4/3)*(xp^(2)+x1^(2)*yp^(2)+x1^(2)*(sin(y1))^(2)*zp^(2)));
 
options=ddeset('AbsTol',1e-15,'RelTol',1e-13);
[s ys]=ode15s(@syseq,1:10000, [ w1;wp;x1;xp;y1;yp;z1;zp],options);
 
 
figure(1);
plot(ys(:,1),ys(:,3),'b');
 
figure(2);
plot(ys(:,1),ys(:,5),'r');
 
 
function dydt = syseq(s,y) 
R0=1;
H0=65000;
cv=3*10^(8);
 
 dydt=zeros(8,1);
 dydt=[ y(2) ;                        
       -((2/3)*R0^(2)*(3*H0/(2*cv))^(4/3)*y(1)^(/3))*y(4)^(2)+y(3)^(2)*y(6)^(2)+y(3)^(2)*(sin(y(5)))^(2)*y(8)^(2) ;
       y(4) ;
       -(4/(3*y(1)))*y(2)*y(4)+y(3)*y(6)^(2)-(sin(y(5)))^(2)*y(8)^(2);                            
       y(6) ;
       -(2/y(3))*y(4)*y(6)-(4/(3*y(1)))*y(2)*y(6)+(sin(2*y(5)))*(y(8)^(2)/2);                           
       y(8) ;
       -2*y(8)*((2/(3*y(1)))*y(2)+(1/tan(y(5)))*y(6)+(1/y(3))*y(4)) ];                   
 
 end
 
end

avec d(theta)/ds=yp=0.001, j'obtiens avec le programme en C le graphique geo2.pdf. Avec le code ci-dessus en matlab, j'obtiens le graphique geo3.pdf qui est incorrect (le tracé en bleu sur geo3.pdf et le tracé en pointillé sur geo2.pdf sont censés être identiques).

Toute suggestion est la bienvenue, merci.

[FR119492]

Salut!

quelles solutions me reste t-il pour corriger les erreurs données avec ode45 ?

Essaie avec ode23s.
Jean-Marc Blanc

[fab13]

ode23s produit les mêmes résultats que ode45. Comment préciser dans l'expression

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
 
[s ys]=ode23s(@syseq,1:10000, [ w1;wp;x1;xp;y1;yp;z1;zp],options);

que ys(:,3) (rho dans le systeme d'ode) doit rester positif ?

[FR119492]

Salut!

que ys(:,3) (rho dans le systeme d'ode) doit rester positif ?

Il est vraisemblable que tu n'intègres pas le même système différentiel en C et dans Matlab. Je pense que, dans le sous-programme qui calcule les dérivées des variables d'état, tu dois avoir un if qui aiguille le calcul de manière à ce que la condition ci-dessus soit toujours satisfaite.
Jean-Marc Blanc