Runge Kutta en Fortran Résolu/Fermé

Question

Bonjour,

pour un calcul de concentration d'éléments radioactifs, je dois faire un programme en fortran en utilisant la méthode de Runge-Kutta. Le problème, c'est que le prof ne nous a donné la méthode que pour une équation, alors que je dois résoudre un gros système.
Il commence comme ça (x, y et z étant les concentrations des éléments à l'instant t)

dx/dt = -ax
dy/dt = ax - by
dz/dt = by +cz

Il faout mettre les equations sous forme de matrice, ce qui est facilement faisable, mais comment appliquer Runge-Kutta sur cette matrice ?

Sacabouffe · Accepted Answer

Salut
Tu appliques de la même manière que d'habitude.
Pour RK4 avec un pas de discrétisation h:

  ( x )
v=( y )
  ( z )

v'(t)=f(t,v(t))

f(t,v(t))=A(t)v(t)

        (  -a(t)       0       0   )
A(t) =  (   a(t)     -b(t)     0   )
        (    0        b(t)    c(t) )

v(t0)=v0

v_(n+1) = v_n + h/6*(k1+2*k2+2*k3+k4)

k1 = f(t_n,v_n) = A_n*v_n

k2 = f(t_n+h/2,v_n+h/2*k1) = A_(n+1/2)*(v_n+h/2*k1) = A_(n+1/2)*(v_n+h/2*A_n*v_n)

k3 =  f(t_n+h/2,v_n+h/2*k2) = A_(n+1/2)*(v_n+h/2*k2) = A_(n+1/2)*(v_n+h/2*A_(n+1/2)*(v_n+h/2*A_n*v_n))

k4 = f(t_n+h,v_n+h*k3) =A_(n+1)*(v_n+h*k3) = A_(n+1)*(v_n+h*A_(n+1/2) *(v_n+h/2*A_(n+1/2)*(v_n+h/2*A_n*v_n))) Donc 
v_(n+1) = v_n + h/6*(A_n*v_n+2*A_(n+1/2) *(v_n+h/2*A_n*v_n)+2*A_(n+1/2)*(v_n+h/2*A_(n+1/2) *(v_n+h/2*A_n*v_n))+A_(n+1)*(v_n+h*A_(n+1/2) *(v_n+h/2*A_(n+1/2)*(v_n+h/2*A_n*v_n)))) Cela dit, si tes coefficients a b c sont constants, une résolution numérique est inutile.

Bonne journée

Dedrak · Answer

Bonjour,
tout d'abord merci pour votre aide.

après avoir bien avancé sur mon programme, je me suis rendu compte que les résultats obtenus étaient faux.

Comment retranscrire un A_(n+1/2) dans le programme, alors que par définition A n'est connu qu'à des pas de temps précis, et donc à des n entiers ?

Faut-il considérer que A_(n+1/2)=A_n ?

Merci d'avance pour votre aide.

Sacabouffe · Answer

Salut
De rien ;-)

En fait ce qu'il faut faire, c'est calculer la solution v sur 2 fois moins de pas de temps que pour A (sur un pas de temps 2 fois plus grand si tu préfères), il faut considérer la formule telle quelle est et ne pas faire l'approximation A_(n+1/2)=A_n.

Tu as A_0, A_1, A_2, ..., A_(2n), A_(2n+1), A_(2n+2), A_(2n+3), ...  et tu calcules v_0, v_2, ..., v_(2n), v_(2n+2), ...

Si tu veux mieux comprendre, tu poses B_n=A_(2n) et w_n=v_(2n) et tu appliques les formules du message 1.
B est alors connu aux pas de temps entiers et demi-entiers et tu calcules w aux pas de temps entiers.

Bonne journée

dedrak · Answer

Merci beaucoup. En fait, le problème c'est que les calculs de A à l'instant n sont dépendants de la valeur de v à ce même pas de temps, sinon ce serait trop facile ^^

C'est justement pour ça que j'avais pris l'hypothèse absurde que A ne changeait pas entre n n+1/2.

Mais merci quand même pour l'aide, je vais devoir me tourner vers une autre méthode qui ne posera pas ce problème pour résoudre le système d'équations.

Sacabouffe · Answer

Non, la méthode de Runge-Kutta fait intervenir A aux pas de temps et aux demi-pas de temps et v aux pas de temps.
https://fr.wikipedia.org/wiki/Méthodes_de_Runge-Kutta
Bonne soirée

Dedrak · Answer

Je sas bien, c'est justement pour ça que j'avais précisé que mon hypothèse était absurde.

Cependant, pour résoudre le calcul, je suis obligé de la faire, puisque je n'ai pas de moyen de calculer plus de valeurs de A que de valeurs de v, les deux étant connectés en permanence.

Sacabouffe · Answer

Oui mais t'as dit je vais devoir me tourner vers une autre méthode qui ne posera pas ce problème pour résoudre le système d'équations.
Il y a pas besoin de se tourner vers une autre méthode, la méthode RK4 convient très bien même si A et v sont connectés en permanence, sauf qu'à la différence d'une méthode d'Euler par exemple, le calcul de v_(n+1) passe par le calcul de A faisant intervenir des demi-pas de temps.
Dans ton message initial t'as pas précisé que les coefficients de la matrice dépendait de v, on savait même pas s'ils dépendaient du temps, alors je ai écrit les équations dans le cas intermédiaire (ni trop simple où les coefficients sont constants, ni trop compliqués où les coefficients dépendent de v).
Si les coefficients dépendent de v, nul besoin de faire intervenir v au demi-pas de temps, A_(n+1/2) s'exprime à partir de v_n de la façon dont c'est expliqué sur le lien Wikipédia que j'ai donné précédemment.
http://fr.wikipedia.org/wiki/Méthodes_de_Runge-Kutta
Bonne nuit

Runge Kutta en Fortran

7 réponses

Discussions similaires