Résolution numérique d'équa. diff.

[Coco44]

Bonjour à tous,

Il y a quelque semaine de cela, on nous a donné un examen d'informatique pour la physique avec en gros le sujet suivant :
On a un pendule pesant de masse m pendu à un fil de longueur R.
Ce pendule fait un angle avec la verticale.
Une étude physique du problème nous donne l'équation différentielle :
... avec :

On ne sais pas résoudre cette équation pour trouver , donc on va utiliser une méthode incrémentale (pas à pas) pour calculer et ses dérivées ( et ).

On nous donne alors :



(formules obtenues par développement de Taylor à l'ordre 2)

C'est là que les soucis commencent :
Le jour du TP, j'ai réussi à avoir de bonnes valeurs : oscillait gentiment entre deux valeurs constantes, avec la bonne fréquence.
Aujourd'hui, en voulant refaire la même chose chez moi, j'obtiens presque la même chose, à ceci près que "diverge" : l'enveloppe de la courbe est une fonction affine (~0.0021*t + ) au lieu d'être constante...

(bleu : résultat obtenu, rouge : résultat "théorique", jaune : enveloppe)
J'ai beau tourner la chose dans tous les sens : je ne vois pas ce que j'ai pu faire pour avoir un résultat si différent.

Voici le code que j'utilise pour calculer les valeurs :

float fT = 0, fY1 = 0.017f, fY2 = 0, fY3 = (-9.81f/0.25f)*sin(0.017f);
// theta à t=0 vaut ~1 degré
for (uint ui = 0; ui < uiPointNbr; ++ui)
{
    // le temps
    fT  += fDT;
    // theta
    fY1 += fDT*fY2 + (fDT*fDT/2)*fY3;
    // theta prime
    fY2 += fDT*fY3;
    // theta seconde
    fY3 = (-9.81f/0.25f)*sin(fY1);
}

En utilisant une autre équation différentielle plus simple :
... j'obtiens un résultat bien plus fidèle à la réalité, mais une erreur subsiste (polynôme du second degré en t).
Pourtant, le développement de Taylor est exact ici, puisque est nul. Il y a donc une erreur dans mon traitement, mais je n'arrive pas à voir où...

Est-ce que vous auriez une piste ?
Merci d'avance !
-----

[KerLannais]

Salut,

J'ai testé ton code sous scilab, il est effectivement instable dans le cas du pendule mais en augmentant le nombre de points l'amplification s'estompe ce qui veux dire que le code est tout à fait correct. De plus j'ai fais aussi le test avec

et j'obtiens la solution exacte à la précision machine près, et donc il n'ya aucun problème avec ton code, peut-être t'es tu trompé quand tu as calculé la solution exacte à la main dans le cas de la dérivée seconde constante.

[Coco44]

Hum, je me serais attendu à ce que le résultat soit nettement plus précis...
J'ai utilisé un incrément de temps égal à 0.01s, c'est peut être trop pour une étude correcte ?

PS : Au sujet de cas simple avec y''' = cte, au temps pour moi : j'avais fais un arrondi dans mes calculs "théoriques" Ça confirme ce que tu disais.

[KerLannais]

Re,

Oui c'est le problème de la stabilité des schémas numérique, c'est la raison pour laquelle en analyse numérique on utilise des schémas plus perfectionnés. En général, les problèmes de stabilité se gèrent plus ou moins en implicitant partiellement ou complètement le schéma.

Je me propose de t'expliquer avec les mains ce qu'on appelle la stabilité d'un schéma.

En dimension 1 tu travail en général sur un segment [a,b] subdivisé, c-à-d qu'on se donne

notons la solution exacte, le schéma te donne un moyen de calculer des approximations de la solution exacte au points de la subdivision. Autrement dit:

En général également le schéma est donné par une récurrence d'ordre 1, c-à-d qu'on a une fonction telle que:

et on pose

(Rq: dans ton exemple il faut prendre pour un couple de valeurs l'une étant l'approximation de et l'autre de , la fonction va donc de dans lui même)

on appelle erreur de convergence la quantité:

En général on décompose cette erreur en deux erreurs: l'erreur de cohérence et l'erreur de stabilité.
L'erreur de cohérence est due au fait qu'à chaque étape on commet une erreur quand on calcule la valeur en en fonction de la valeur en , en effet, quand bien même on supposerait que, au miracle, la valeur est exactement égale à la valeur exacte on aurait quand même une erreur sur

Puisque qu'en générale la fonction est une approximation discrète du problème continu dont est solution et qu'en plus il peut y avoir des erreurs d'arrondis. On définit donc l'erreur de cohérence par:

En général la fonction dépend du pas de la subdivision et plus ce pas est petit, plus approche fidèlement le problème continu. Si l'erreur de cohérence tend vers 0 quand le pas tends vers 0 on dit que le schéma est cohérent.
On a:

Autrement dit, l'erreur de convergence au point est plus petite que l'erreur de cohérence plus une erreur qu'on appelle erreur de stabilité:

c'est en quelque sorte l'image par de l'erreur de convergence au point :

En fait, si

on dit que le schéma est stable (Rq: il y a plusieurs définition de stabilité dans la litérature). Cela siginfie que le schéma n'amplifie pas à chaque étape les erreurs déjà commises (qui sont le résultat du cumul des erreurs de cohérence à chaque étape). Ainsi, quand un schéma est stable, pour que le schéma converge il suffit que l'erreur de cohérence à chaque étape soit négligeable devant le pas ( c-à-d ) quand le pas tends vers 0 de sorte que la somme des erreurs de cohérence tende vers 0 quand le pas tends vers 0. C'est le fameux principe "convergence=stabilité+cohéren ce".

[A voir en vidéo sur Futura]

[Coco44]

J'ai tout saisi, sauf :

En fait, si

on dit que le schéma est stable...

... c'est pas plutôt ? (ce qui reviendrait à dire que l'erreur de stabilité est décroissante ou constante).
En tout cas c'est une très bonne explication, merci !

En calculant la valeur de ces erreurs dans le cas du pendule, j'en arrive à la conclusion que le schéma est cohérent, mais pas stable.

[KerLannais]

Oui c'est ça

désolé pour la faute de frappe.