résolution numérique équation différentielle

[Denizou]

Bonjour,

je souhaiterais résoudre numériquement une équation différentielle y'=f(x,y) avec les méthodes d'Euler et du point milieu en xf avec (x0,y0) les conditions initiales.
Dans un 1er temps, j'ai fixé le nombre d'opération à n, et j'ai écrit une boucle for(k:=1;k<=n;k++) suivi des calculs avec h=(xf-x0)/n. Pas de problème.

Dans un 2nd temps, j'ai voulu fixé h et faire varier xf pour obtenir un graphique.
1er essai : j'ai écrit une boucle while(x!=xf) en faisant attention à ce que le pas permette d'obtenir xf. Lorsque h est très petit (0,001 par exemple) il n'y a pas d'arrêt
2ème essai : j'ai écrit une boucle while(x<=xf). Le programme s'arrête mais les valeurs ne sont pas correctes selon les valeurs de xf et plus bizarre, le point milieu converge moins vite qu'Euler !
3ème essai : je calcule n=(xf-x0)/h et j'écris une boucle for. Tout fonctionne. Mais lorsque j'utilise seq pour obtenir ma liste de points il y a un décalage des valeurs...

Quel est le problème rencontré ? S'agit-il d'un problème avec les flottants ?

Pour être plus précis je cherchais à résoudre y'=sin(xy) avec(0,3) comme condition initiale.

Merci de votre aide

[parisse]

L'absence d'arret pour while(x!=xf) me parait plausible, a cause des erreurs d'arrondis. Pour le reste, il faudrait avoir le programme sous les yeux, difficile de donner des elements de reponse a priori.

[Denizou]

Ce qui m'interpelle le plus :

Code : Tout sélectionner

Euldif(f,x,y,xf,h):={
  local k,n;
  n:=floor((xf-x)/h);
 for(k:=1;k<=n;k++){
    y:=y+h*f(x,y);
    x:=x+h;  }
    retourne evalf(y,6);
}:;

Ensuite

Code : Tout sélectionner

f(x,y):=sin(x,y)

Puis

Code : Tout sélectionner

seq(Euldif(f,0,3,k,0.1),k=0..2,0.1)

donne pour résultat

Code : Tout sélectionner

(3.0,3.0,3.02955,3.02955,3.0865,3.16642,3.26183,3.36165,3.45186,3.51819,3.55032,3.54495,3.5057,3.44016,3.35675,3.26276,3.16381,3.06386,2.96566,2.87103,2.78122)

Euldif(f,0,3,0.3,0.1) renvoie bien 3.0865 au lieu de 3.02955 affiché plus haut. Il y a donc un décalage des valeurs. Pourquoi ?

Pour la boucle avec while :

Code : Tout sélectionner

Euler_dif(f,x,y,xf,h):={
     while(x<xf){
    y:=y+h*f(x,y);
    x:=x+h;}
  
  retourne evalf(y,6)
}:;

En fait en regardant à tête reposée et avec le debug, l'inégalité stricte n'est pas prise en compte c'est ce qui pose problème.

Merci de votre aide.

[parisse]

Je viens d'essayer sur mac, je n'observe pas de decalage. Mais ca depend peut-etre de l'architecture. Je pense qu'il serait interessant de renvoyer non pas y mais [x,y] a la fin de la boucle, ca permet de savoir en quel abscisse on a une valeur approchee de y (meme chose pour le 2eme programme avec while).

[Denizou]

Je viens d'essayer sur mac, je n'observe pas de decalage. Mais ca depend peut-etre de l'architecture. Je pense qu'il serait interessant de renvoyer non pas y mais [x,y] a la fin de la boucle, ca permet de savoir en quel abscisse on a une valeur approchee de y (meme chose pour le 2eme programme avec while).

Alors avec la boucle for j'obtiens 2 fois la valeur 0,2
[0.1,3.0],[0.2,3.02955],[0.2,3.02955],[0.3,3.0865]

Avec la boucle while j'obtiens bien les bons résultats mais si a:=odesolve(sin(x*y),[x,y],[0,3],2)
Euler_dif(f,0,3,2,0.001)[1]-a renvoie -0.00038167556681
Milieu_dif(f,0,3,2,0.001)[1]-a renvoie -0.000827761605706

par contre
Euler_dif(f,0,3,2,2/10204.)[1]-a renvoie -7.48801463715e-05
Milieu_dif(f,0,3,2,2/1024.)[1]-a renvoie -7.26065877643e-07
ce qui est plutôt rassurant !

Qu'est-ce qui explique cela ?

Merci de votre aide

[parisse]

C'est surement la valeur finale de x qui explique le 1er cas. Comme 0.001 est represente en base 2, il y a arrondi/troncature et si vous affichez la valeur finale de x avec Euler_dif(f,0,3,2,0.001)[ vous verrez que ce n'est pas 2, mais 2.001 avec boucle while. Le probleme n'apparait plus pour un pas represente exactement en base 2 comme 2/1024.
Donc mon conseil c'est d'utiliser une boucle for, d'ajouter un 1e-10 dans le floor pour les erreurs d'arrondi et de reajuster le pas pour que xf=x+n*h
n:=floor((xf-x)/h+1e-10);
h:=(xf-x)/n;

[Denizou]

C'est surement la valeur finale de x qui explique le 1er cas. Comme 0.001 est represente en base 2, il y a arrondi/troncature

Je me doutais bien qu'il devait y avoir un problème de ce genre...

Donc mon conseil c'est d'utiliser une boucle for, d'ajouter un 1e-10 dans le floor pour les erreurs d'arrondi et de reajuster le pas pour que xf=x+n*h
n:=floor((xf-x)/h+1e-10);
h:=(xf-x)/n;

Mais pour la boucle while le problème vient que l'inégalité stricte n'est pas respectée...
D'ailleurs que je mette < ou <=, j'ai exactement les mêmes résultats tant pour x que pour y.
J'ai utilisé le debug avec un pas de 0.1 pour le calcul en xf=0.5. Lorsque Xcas arrive à 0.5, il ne s'arrête pas et fait une boucle de plus : il renvoie donc la valeur en xf=0.6

Ce que je ne comprends pas c'est qu'en utilisant seq il calcule correctement....

Du coup, pour obtenir le bon résultat le test d'arrêt doit être x<xf-h mais c'est tiré par les cheveux. De plus, avec ce test je ne pourrai pas travailler avec seq pour avoir des résultats cohérents puisqu'alors j'aurais deux fois la valeur en 0 et ma liste s'arrêtera à 1.9 au lieu de 2.

Merci de votre aide.

[parisse]

Utiliser un while avec des données approximatives c'est dangereux, à cause des erreurs d'arrondi, on risque en effet de faire une itération de plus. C'est pour cela qu'il vaut mieux faire un for comme je l'ai proposé. Ensuite plutot que de faire un seq d'appels à Euler, il est beaucoup plus efficace de stocker dans une liste les valeurs de x et de y intermédiaires, comme ça vous n'appelez votre programme qu'une seule fois, et ces valeurs ne sont calculées qu'une seule fois.

P.S.: qui n'a rien à voir avec le fond, pouvez-vous utiliser l'adresse de forum xcas.e.ujf-grenoble.fr au lieu de pcm1.e.ujf-grenoble.fr?

[Denizou]

Utiliser un while avec des données approximatives c'est dangereux, à cause des erreurs d'arrondi, on risque en effet de faire une itération de plus.

Je ne pensais pas que cela poserait problème avec un pas de 0.1.

Ensuite plutot que de faire un seq d'appels à Euler, il est beaucoup plus efficace de stocker dans une liste les valeurs de x et de y intermédiaires, comme ça vous n'appelez votre programme qu'une seule fois, et ces valeurs ne sont calculées qu'une seule fois.

Merci pour cette astuce et pour vos réponses toujours aussi rapides.
Pour info, les algorithmes de résolution d'équation différentielles sont proposés avec une boucle while et l'arrêt à b-h sur le site "Algorithmique et traduction pour Xcas" de Renée De Graeve

P.S.: qui n'a rien à voir avec le fond, pouvez-vous utiliser l'adresse de forum xcas.e.ujf-grenoble.fr au lieu de pcm1.e.ujf-grenoble.fr?

Normalement les changements sont pris en compte.

[parisse]

Je ne pensais pas que cela poserait problème avec un pas de 0.1.

Parce que vous pensez en base 10, comme tout le monde ... sauf les ordis! Il ne faut pas oublier que les seuls rationnels représentables exactement en base 2 ont un dénominateur qui divise une puissance de 2 (et non 10).

Merci pour cette astuce et pour vos réponses toujours aussi rapides.
Pour info, les algorithmes de résolution d'équation différentielles sont proposés avec une boucle while et l'arrêt à b-h sur le site "Algorithmique et traduction pour Xcas" de Renée De Graeve

Merci de le signaler, j'en discuterai avec Renée!
a+

[Denizou]

Juste pour terminer cette discussion, le plus simple, me semble - t - il, c'est encore de demander n puis de calculer le pas. Du coup, on voit bien les résultats obtenus pour x selon que l'on choisit une puissance de 2 ou de 10 et tous les résultats sont cohérents !