Discussion :

[Suzu67]

Bonjour à tous,

Je dois adapter le Schéma d'Euler à l'équation suivante: f(t,u)=(1-t)/(1+u)

Voici mon code Python mais j'obtiens des valeurs bizarres à partir de 2,5.

Quelqu'un pourrait-il m'expliquer où est mon erreur?

On m'a demander dans une précédente question de faire une représentation graphique de U(t) solution de l'équation différentielle y'=(1-t)/(1+y) sur l'intervalle [0;3] donc j'ai pris le même intervalle mais peut-être aurais-je dû en prendre un autre?

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
def fun(t,u):
    return (1-t)/(1+u)
# parametres
t0,tf,y0,n=0,3,0,100
temps=np.linspace(t0,tf,n+1)
 
#euler 
def euler(fct,ti,tf,nt,yi):
    y=np.zeros(nt+1)
    h=(tf-ti)/nt
    t=np.linspace(ti,tf,nt+1)
    y[0]=yi
    for i in np.arange(nt):
        y[i+1]=y[i]+h*fct(t[i],y[i])
    return t,y
# appel de la fonction
 
t,y=euler(fun,t0,tf,n,y0) 
plt.plot(t,y)
plt.title ('intégration Euler Explicite')
plt.xlabel('t')

Merci d'avance de votre aide

[lg_53]

Bonjour

Vous semblez avoir déjà pas mal avancé. Je vais vous poser d'autres questions/remarques qui devrait vous amener sur la bonne piste.

Vous avez précédemment calculer la solution analytique :
- Ok donc tracer là sur le même graphique pour comparer les courbes ! Ou à la rigueur avec 1 subplot, à gauche la courbe d'Euler, à droite la courbe exacte. Et postez en le code, comme ca vous montrez en même temps ce qu'on est censé obtenir. Vous ne remarquez rien sur le tracé de cette courbe entre 0 et 3 ? Que vaut cette fonction en 2.5 ? en 3 ? Et python lui il en dit quoi quand à ces valeurs ? Qu'a-t-il calculé ?
- Quelle est l'ensemble de définition de cette solution analytique ?

PS : Pour que l'on puisse exécuter rapidement le code que vous postez, prenez l'habitude de poster également les import nécessaire.

[Suzu67]

Bonjour lg_53 ,

Oui j'ai une courbe que j'ai faite précédemment.

En faisant print de ma solution U(t) pour voir les valeurs et à partir d'un certain moment j'obtiens des valeurs nan.

Je cherche donc à savoir quelle est l'ensemble exacte des solutions de mon équation.

Or je n'arrive à faire un code qui me donne ma limite à droite (à force de test j'ai trouvé 2.414213 que j'ai remis dans mon code comme 2.4).
Et là je n'ai plus de valeurs bizarre sur mon graphique.

Mais n'y a-t-il pas un moyen plus précis pour faire cela (je débute en Python et m'excuse si cette question est évidente)?

Voici mon code avec les 2 courbes à la fin normalement:

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
import numpy as np
import matplotlib.pyplot as plt
 
a=-1
b=2
c=1
d=-1
 
def U(x):
    return np.sqrt(a*x**2+b*x+c)+d
 
def fun(t,u):
    return (1-t)/(1+u)
# parametres
t0,tf,y0,n=0,2.4,0,100
temps=np.linspace(t0,tf,n+1)
 
#euler 
def euler(fct,ti,tf,nt,yi):
    y=np.zeros(nt+1)
    h=(tf-ti)/nt
    t=np.linspace(ti,tf,nt+1)
    y[0]=yi
    for i in np.arange(nt):
        y[i+1]=y[i]+h*fct(t[i],y[i])
    return t,y
# appel de la fonction
 
t,y1=euler(fun,t0,tf,n,y0) 
plt.plot(t,y1)
plt.title ('Intégration Euler Explicite de f(t,u) ')
plt.xlabel('t')
 
 
x=np.linspace(0,2.4)
y=U(x)
plt.plot(x,y, label='U(t)')
plt.plot(t,y1, label='Euler')
plt.title ('Comparaison')
plt.legend()
plt.xlabel('t')

Merci d'avance de votre réponse

[lg_53]

1) L'ensemble de définition se calcule précisément depuis la solution analytique puisque vous l'avez ! Et correspond donc au x tel que -x²+2x+1 >= 0 (étude de signe d'un polynôme du second degré), et ça vous donnera l'ensemble de définition qui devrait être [1-sqrt(2); 1+sqrt(2)]
2) Numériquement, lorsqu'on n'a pas le luxe d'avoir une solution analytique, et bien on analyse la fonction f. Et en l’occurrence, on voit que si u vaut -1, alors f n'est pas définie. Numériquement il y a très peu de chance que l'on tombe pile sur la valeur u=-1, donc les calculs seront possibles. Mais en pratique, la fonction f admettant une discontinuité en u=-1, cela va créé de forte instabilité dans le schéma d'Euler, qui n'est passé au travers de cette discontinuité que parce que la précision numérique fait qu'on n'est pas passé pile par la valeur u=-1, d'où l'apparition d'oscillation.

Il existe des schémas à pas adaptatifs (on peut faire du Euler à pas adaptatif), qui eux ne créerait pas d'instabilité. Il diminuerait de plus en plus leur pas de temps à l'approche de cette discontinuité, car la dérivée deviendrait de plus en plus grande. Jusqu'à ce que le pas de temps en soit réduit au epsilon machine, et là :
- soit votre schéma contient un garde fou pour soi détecter ces cas de figures soi limiter le nombre d'itération de sorte à arrêter la boucle et renvoyer un warning à l'utilisateur
- sinon, votre schéma tournera en boucle, sans jamais s'arrêter


Il peut être intéressant aussi de comparer et d'examiner comment se comporte scipy.ode.integrate sur ce cas précis.