Discussion :

[isaz24]

Bonjour,
je dois utiliser la méthode de runge-kutta pour calculer la pression dans un film de fluide

au final , j'ai une équation qui me donne la dérivée seconde de la pression en fonction de la position et de la dérivée première de la pression :

d²p/dx²=f(x, dp/dx)

j'ai utilisé la méthode de runge kutta d'ordre 4 pour faire le calcul de dp/dx sur toute la longueur de la zone à étudier, mais je ne sais pas comment calculer la pression avec une bonne précision.
j'obtiens bien une assez grande précision sur dp/dx, mais je n'arrive pas à obtenir une bonne précision pour p

donc ma question : comment utiliser les facteurs k1 à k4 de la méthode pour calculer p avec une bonne précision ?

(pour le moment , j'ai commencé par calculer p en utilisant dp/dx, puis j'ai ajouté la dérivée seconde pour passer à un calcul d'ordre 2 , mais maintenant j'aimerais savoir comment prendre en compte les facteurs k pour passer aux ordres supérieurs)

[marco056]

Bonjour,
Sans code, tu auras peu de réponses.

[tyrtamos]

Bonjour,

Juste une idée.

Il peut arriver que dans un calcul mathématique, on ait besoin d'une plus grande précision dans les calculs intermédiaires pour obtenir un résultat final acceptable. Ça m'est déjà arrivé par exemple dans les calculs des nombres de Bernouilli.

Si c'est le cas, on peut "pousser" la précision des calculs intermédiaires en utilisant le module "decimal" (on peut fixer le nombre de chiffres significatifs à 1000 si nécessaire!). S'il y a en plus un problème de puissance de calcul, il y a d'autres solutions plus "musclées" comme avec les modules externes "gmpy2" ou "mpmath". Ces modules externes s'installent facilement avec pip. Mais, bien sûr, dans ces 3 cas, il faut réorganiser le code de calcul.

[isaz24]

pour le moment , mon code ressemble à ça :

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
 
p=0
dpdx=0
 
i=-1
while (i<N-1) :
	i=i+1
	rk1=dx*fonction(x0+i*dx , dpdx) 
	rk2=dx*fonction(x0+(i+0.5)*dx , dpdx+rk1/2)
	rk3=dx*fonction(x0+(i+0.5)*dx , dpdx+rk2/2)
	rk4=dx*fonction(x0+(i+1)*dx , dpdx+rk3)
 
	a=(rk1+rk2*2+rk3*2+rk4)/6
 
	p=p+dx*(dpdx+a/2)
	dpdx=dpdx+a

fonction() calcule la dérivée seconde de p par rapport à x, en fonction de x et de dpdx
donc, j'utilise la méthode de runge-kutta d'ordre 4 pour calculer dpdx sur mon intervalle, et j'obtiens une bonne précision avec peu de points, et une erreur en dx**5.
par contre , pour le calcul de p, je n'arrive pas à obtenir une bonne précision à moint d'augmenter radicalement le nombre de points.
et au mieux, j'ai une précision en dx**1
je pense qu'il doit être possible de calculer p avec une meilleure précision en utilisant les facteurs de runge-kutta plutot que de faire simplement la méthode des trapèzes pour obtenir un résultat avec une précision en dx**5 comme pour dpdx
(pour le moment j'utilise simplement la moyenne de dpdx(i) et dpdx(i+1) pour calculer p à chaque pas)

[Ehouarn]

Bonjour,

fonction() calcule la dérivée seconde de p par rapport à x, en fonction de x et de dpdx

Je ne vois pas comment, vu qu'à ce stade p n'est pas connu, à moins de faire d'énormes approximations (probablement peu justifiables).

Je pense qu'il y a là surtout un souci de méthode. Runge-Kutta est habituellement utilisé pour intégrer dune EDO (d'ordre 1) à partir d'une condition initiale, et tu veux l'utiliser pour résoudre une EDP d'ordre 2 (avec conditions aux limites; tu n'en parles pas, mais j'imagine que c'est le cas).

[isaz24]

Bonjour,


Je ne vois pas comment, vu qu'à ce stade p n'est pas connu, à moins de faire d'énormes approximations (probablement peu justifiables).

Je pense qu'il y a là surtout un souci de méthode. Runge-Kutta est habituellement utilisé pour intégrer dune EDO (d'ordre 1) à partir d'une condition initiale, et tu veux l'utiliser pour résoudre une EDP d'ordre 2 (avec conditions aux limites; tu n'en parles pas, mais j'imagine que c'est le cas).

C'est bien une équation différentielle du premier degré, car j'ai une équation de la forme
dy/dx =f(x,y)

(ici , y=dp/dx )

le seul problème c'est qu'une fois que j'ai résolu l'équation et trouvé y(x) tout le long de mon intervalle, j'aimerais calculer sa primitive Y(x)
donc ici, Y(x)=p

l'équation qui me donne dy/dx =f(x,y) vient des équations de navier-strokes en mécanique des fluides , combinées à l'équation de conservation de la masse dans le système (donc le débit est le même partout),
et en prenant en compte les conditions aux limites dans 2 directions du système.

à la fin, ça me donne :

d²y/dx² = (6*vitesse*mu/h^3 - 3/h*dp/dx)*dh/dx

avec h l'épaisseur du film de fluide, mu la viscosité, et vitesse est la vitesse de déplacement d'une paroi (le film de fluide est entre 2 parois, l'une fixe, et l'autre mobile)

vitesse et mu sont des constantes, h et dh/dx dépendent de x ; et il reste dp/dx

donc on a bien une fonction qui dépend de x et de dp/dx
p n'intervient nulle part dans l'équation, mais je voudrais le calculer car j'en ai besoin pour d'autres calculs
(d'ailleurs j'aurais aussi besoin de calculer l'intégrale de p plus tard , mais c'est moins important et j'aurais moins besoin de précision)

[Ehouarn]

Bonjour,

De ce que je comprends de ton problème, il s'agit bien de résoudre une équation différentielle avec conditions aux limites, or une méthode d'intégration type Runge Kutta part d'un point pour "avancer" (habituellement dans le temps) et tu voudrais l'appliquer en partant d'une limite du domaine jusqu'à l'autre borne. Mais en tentant de faire cela tu n'utilises aucune information sur la condition à la limite sur la borne opposée. Donc, à moins d'un miracle, ta solution ne sera pas "raccord" avec cette condition, ce qui illustre bien inapplicabilité de la méthode.

Pour résoudre un problème avec conditions aux limites, on discrétise plutôt l'équation à résoudre (par exemple avec une méthode de différences finies) qui permet de relier les valeurs aux points de maillage entre eux et leur relation avec les valeurs aux limites. Ce qui amène alors à écrire un système matriciel à résoudre ensuite par une méthode numérique adaptée.

[isaz24]

en fait, j'arrive assez bien à résoudre le problème.
pour prendre en compte les conditions aux limites aux deux bouts (pression nulle aux deux bouts), je fais un ajustement :
je pars de la gauche de mon domaine avec la pression = 0, et une valeur arbitraire pour sa dérivée, j'utilise runge-kutta pour calculer la pression à droite.
j'obtiens une valeur non-nulle, donc différente de ma condition limite,
puis je refais le calcul en repartant de la gauche avec pression(gauche) = 0, et une autre valeur initiale pour la dérivée de la pression (à choisir .

mon ajustement arrive en général à trouver un bon résultat après à peu près 4 essais (je tombe sur une pression inférieure à 10^-6 à droite , alors que les pressions dans le système sont de l'ordre de 10^4 )


mais c'est vrai que la méthode des différences finies a l'air intéressante, je crois que je vais essayer d'utiliser ça.


le seul problème , c'est que je suis obligé de découper l'espace en un grand nombre de points pour obtenir une bonne précision sur la pression, alors que j'arrive à atteindre une très bonne précision sur dpdx avec bien moins de points.

[Ehouarn]

Bonjour,

en fait, j'arrive assez bien à résoudre le problème.
pour prendre en compte les conditions aux limites aux deux bouts (pression nulle aux deux bouts), je fais un ajustement :
je pars de la gauche de mon domaine avec la pression = 0, et une valeur arbitraire pour sa dérivée, j'utilise runge-kutta pour calculer la pression à droite.
j'obtiens une valeur non-nulle, donc différente de ma condition limite,
puis je refais le calcul en repartant de la gauche avec pression(gauche) = 0, et une autre valeur initiale pour la dérivée de la pression (à choisir .

mon ajustement arrive en général à trouver un bon résultat après à peu près 4 essais (je tombe sur une pression inférieure à 10^-6 à droite , alors que les pressions dans le système sont de l'ordre de 10^4 )

Je dirai que tu as une sacré veine que cette magouille fonctionne.
[EDIT]
Je ne devrai pas être aussi expéditif, c'est en fait la méthode du tir que tu utilises là; je l'avais oublié (Alzeihmer, déjà?).
[/EDIT]

le seul problème , c'est que je suis obligé de découper l'espace en un grand nombre de points pour obtenir une bonne précision sur la pression, alors que j'arrive à atteindre une très bonne précision sur dpdx avec bien moins de points.

Là je suis à nouveau un peu perdu. Si tu disposes déjà par ailleurs des valeurs de dpdx en un certain nombre de points alors il ne s'agit plus de résoudre une équation différentielle mais seulement de faire une intégration numérique avec une méthode du type trapèzes (ou plus évoluée si la précision est importante).

[isaz24]

oui, c'est vrai, c'est plutot ça qu'il me faudrait en fait. j'aimerais pouvoir calculer l'intégrale avec une méthode semblable à celle des trapèzes mais en plus évoluée.
a mimimum, avec la méthode des trapèzes, je peux utiliser la dérivée première que je connais, pour avoir une méthode d'ordre 1
p(i+1)=p+dx*dpdx

comme j'ai une fonction qui me donne la dérivée seconde, je peux l'utiliser pour obtenir une méthode d'ordre 2
p(i+1)=p+dx*dpdx+dx*dx/2*f(x,dpdx)

comme dans la méthode de runge kutta, on a k1=dx*f(x,dpdx), je remplace :

p(i+1)=p(i)+dx*dpdx+dx/2*k1


et j'espérais sur le même principe pouvoir utiliser les autres coefficients k pour passer aux ordres supérieurs.

[Ehouarn]

Bonjour,

Formellement une méthode de RKn est d'ordre n parce que le schémas coïncide avec le développement de Taylor à l'ordre n de la fonction recherchée. Donc oui, en réécrivant ce développement tu devrais pouvoir retrouver les coefficients et donc intégrer "proprement" la fonction à l'ordre n.

Une autre approche pourrait être d'utiliser une méthode d'intégration du type Newton-Cotes (une généralisation de la méthode des trapèzes) d'ordre approchant celui du RK, par exemple la formule de "Boole-Villarceau" (exacte pour un polynôme de degré <=5). L'avantage est alors de ne se baser que sur les points successifs obtenus par l'intégration via RK, sans se préoccuper des intermédiaires.