Transcript Document 4597942

Informatique en PCSI et MPSI
Champollion 2013-2014
Méthodes d’Analyse Numérique
Implémentation et Application en Python
Équations différentielles ordinaires
A. HASSAN
19 février 2014
A. Hassan@Champollion
PCSI-MPSI – 1
Résolution des équation différentielles
ordinaires (EDO)
Résolution des équation
différentielles ordinaires
(EDO)
Le Problème
Le Problème (suite)
Problème de Cauchy
Méthodes de résolution:
Ordres 1 et 2
Méthode d’Euler ou RK1
(Runge-Kutta d’ordre 1)
Implémentation de
l’algorithme d’Euler
Implémentation de la
méthode d’Euler en
Python
Runge Kutta d’ordre 4
(RK4)
Méthode de Runge-Kutta
d’ordre 4 en Python
Méthode de Runge-Kutta
d’ordre 4 en Python
Utilisation de la
commande odeint du
module
scipy.integrate
Commande
odeint (suite)
A. Hassan@Champollion
Exercice: L’équation de
Van Der Pol (1924)
Implémentation en Python
Exploitation graphique
Exploitation graphique:
PCSI-MPSI – 2
Champs de vecteurs
Solution
Le Problème
Soit I un intervalle de R, D un ouvert de R n et f une fonction connue de
R n × I −→ R n . On cherche y : I → R n telle que




y1 ( t )
f 1 ( y1 ( t ), · · · , y n ( t ); t )
 . 


..
et
y
(
t
)
=
y′ (t) = f (y(t); t) où f (y(t); t) = 
 .. 

.
f n ( y1 ( t ), · · · , y n ( t ); t )
yn (t)
Ce type d’équations s’appelle Équation résolue en y′ .
Par contre sin(t.y′ + cos(y′ + y)) = t.y + y′ n’est pas résolue en y′ .
I.e. Impossible d’écrire y′ = f (t, y) (avoir y′ explicitement en fonction de t et y)
Ordinaire : la dérivation est effectuée par rapport à t uniquement.
∂2 T
∂2 T
∂T
= 2 + 2 , fonction inconnue
Équation aux dérivées partielles(EDP) :
∂t
∂x
∂y
T (t, x, y) (Équation de la chaleur ou de la diffusion)
A. Hassan@Champollion
PCSI-MPSI – 3
Le Problème (suite)
Problème de Cauchy (ou dela condition initiale) :
y1 ( t )


Trouver t 7→ y(t) =  ...  de I −→ R n telle que
yn (t)
y ( t0 ) = y0
y ′ ( t ) = f ( y ( t ); t )
Problème des conditions aux limites :
Cas où l’on dispose d’informations sur y à des valeurs différentes de t.
Conditions initiales
 ′
 y1 (t) = −y1 (t) + y2 (t) + sin(t)
y2′ (t) = 2y1 (t) − 3y2 (t)

(y1 (0); y2 (0)) = (0, 3)
A. Hassan@Champollion
Conditions
aux limites
 ′
y (t) = −y1 (t) + y2 (t) + sin(t)


 1′
y2 (t) = 2y1 (t) − 3y2 (t)
y1 (0) = 2



y2 ( π4 ) = −1
PCSI-MPSI – 4
Problème de Cauchy
Sauf quelques cas particuliers, il est presque impossible de
trouver des solutions analytiques(i.e. se ramener à un calcul de
primitives).
On cherche des solutions approchées, donc Résolutions
numériques.
Principe : Si l’on connaît y à l’instant (abscisse) t, comment
obtenir y à l’instant t + h, y(t + h), puis recommencer avec
y(t + h) ?
Choisir h assez petit et utiliser les Développements limités
y′ ( t )= f ( t,y(t)
1 2 ′′
n
y(t + h) = y(t) + hy (t) + h y (t) + · · · + o(h ) ===========⇒
2
′
= y(t) + Φ(t, h, f (y, t)) + o(hn )
On cherche Φ(t, h, f (y, t)) de telle sorte l’ordre de o(hn ) soit le
plus élevé possible (en minimisant le coût et la complexité des
calculs).
A. Hassan@Champollion
Résolution des équation
différentielles ordinaires
(EDO)
Le Problème
Le Problème (suite)
Problème de Cauchy
Méthodes de résolution:
Ordres 1 et 2
Méthode d’Euler ou RK1
(Runge-Kutta d’ordre 1)
Implémentation de
l’algorithme d’Euler
Implémentation de la
méthode d’Euler en
Python
Runge Kutta d’ordre 4
(RK4)
Méthode de Runge-Kutta
d’ordre 4 en Python
Méthode de Runge-Kutta
d’ordre 4 en Python
Utilisation de la
commande odeint du
module
scipy.integrate
Commande
odeint (suite)
Exercice: L’équation de
Van Der Pol (1924)
Implémentation en Python
Exploitation graphique
Exploitation graphique:
PCSI-MPSI – 5
Champs de vecteurs
Solution
Méthodes de résolution : Ordres 1 et 2
• Le plus souvent h = (t f − t0 )/N où [t0 ; t f ] désigne
l’intervalle de résolution et N le nombre de sous-intervalles.
• Sinon : On considère T = {t0 ; t1 ; · · · ; tn−1 ; tn } et
h i = ( t i +1 − t i )
Ainsi y(t) est l’estimation courante en t et y(t + h) l’estimation
au pas suivant
Pour n = 1 ⇒
Méthode d’Euler (ou Runge Kutta d’ordre 1).
y(t + h) ≈ y(t) + h f (t, y(t))
Pour n = 2 ⇒
Méthode de Runge Kutta d’ordre 2.
K1 = h. f (t, y(t))
1 1
K2 = h. f (t + h; K1 + y(t))
2 2
y ( t + h ) ≈ y ( t ) + K2
A. Hassan@Champollion
Résolution des équation
différentielles ordinaires
(EDO)
Le Problème
Le Problème (suite)
Problème de Cauchy
Méthodes de résolution:
Ordres 1 et 2
Méthode d’Euler ou RK1
(Runge-Kutta d’ordre 1)
Implémentation de
l’algorithme d’Euler
Implémentation de la
méthode d’Euler en
Python
Runge Kutta d’ordre 4
(RK4)
Méthode de Runge-Kutta
d’ordre 4 en Python
Méthode de Runge-Kutta
d’ordre 4 en Python
Utilisation de la
commande odeint du
module
scipy.integrate
Commande
odeint (suite)
Exercice: L’équation de
Van Der Pol (1924)
Implémentation en Python
Exploitation graphique
Exploitation graphique:
PCSI-MPSI – 6
Champs de vecteurs
Solution
Méthode d’Euler ou RK1 (Runge-Kutta d’ordre 1)
Principe de la méthode :
Chercher une solution approchée de
t f − t0
y′ = f (t; y)
sur l’intervalle I = [t0 ; t f ]
y0 = y ( t0 )
) où
h = ( t i +1 − t i ) ⇒ y ( t i + h ) ≈ y ( t i ) + h i . f ( y ( t i ), t i )
Géométriquement : Remplacer la courbe sur [ti ; ti + h] par la tangente.
Si Mn (tn ; yn ) est le point courant de la courbe "solution", le nouveau point :
Mn+1 (tn + h; yn + h f (tn , yn )).
pas d’intégration : h = (
N
b
δy
bc
yn
erreur commise
bc
hn
A. Hassan@Champollion
bc
bc
tn
t n +1
PCSI-MPSI – 7
Algorithme d’Euler (Runge Kutta d’ordre un)
Euler ( f , y0 , t0 , t f , N )
Entrées :
f fonction données
(t0 ; y0 ) point initial
t f abscisse final
N nombre de points de [t0 ; t f ]
Sorties : Ly liste des ordonnées yk , k = 0; 1; · · · ; N
( t f − t0 )
h←
N
L y ← y0 , L t ← t0
pour k de 1 à N faire
y0 ← y0 + h. f (t0 , y0 )
t0 ← t0 + h
L y ← L y , y0 ;
L t ← L t , t0
retourner Ly et Lt
A. Hassan@Champollion
# stocker les solutions
# stocker les abscisses
Résolution des équation
différentielles ordinaires
(EDO)
Le Problème
Le Problème (suite)
Problème de Cauchy
Méthodes de résolution:
Ordres 1 et 2
Méthode d’Euler ou RK1
(Runge-Kutta d’ordre 1)
Implémentation de
l’algorithme d’Euler
Implémentation de la
méthode d’Euler en
Python
Runge Kutta d’ordre 4
(RK4)
Méthode de Runge-Kutta
d’ordre 4 en Python
Méthode de Runge-Kutta
d’ordre 4 en Python
Utilisation de la
commande odeint du
module
scipy.integrate
Commande
odeint (suite)
Exercice: L’équation de
Van Der Pol (1924)
Implémentation en Python
Exploitation graphique
Exploitation graphique:
PCSI-MPSI – 8
Champs de vecteurs
Solution
Implémentation de l’algorithme d’Euler
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
1.
On charge les modules numpy (ou scipy ) pour les calculs
2.
3.
numériques et matplotlib.pyplot (ou pylab ) pour les
graphiques.
La fonction (y; t) 7→ f (y; t) est définie au préalable.
On rappelle que y0 = y(t0 ) ∈ R n éventuellement
import numpy as np
import matplotlib . pyplot as mp
##
def Euler (f , t0 , y0 , tf ,N ):
h =.........
# definiton du pas fixe
# sinon h =t [ i +1] - t[ i ]
Ly , Lt =[...] ,[...]
# stockage du point initial
f o r k in ..........: # A la limite Lt peut etre donne ailleurs
# prevoir que y0 peut etre vectoriel
y0 =...................... ...
t0 +=..................... ...
................... ... ... ...
................... ... ... ...
r e t u r n array ( Lt ) , array ( Ly )
##
A. Hassan@Champollion
Résolution des équation
différentielles ordinaires
(EDO)
Le Problème
Le Problème (suite)
Problème de Cauchy
Méthodes de résolution:
Ordres 1 et 2
Méthode d’Euler ou RK1
(Runge-Kutta d’ordre 1)
Implémentation de
l’algorithme d’Euler
Implémentation de la
méthode d’Euler en
Python
Runge Kutta d’ordre 4
(RK4)
Méthode de Runge-Kutta
d’ordre 4 en Python
Méthode de Runge-Kutta
d’ordre 4 en Python
Utilisation de la
commande odeint du
module
scipy.integrate
Commande
odeint (suite)
Exercice: L’équation de
Van Der Pol (1924)
Implémentation en Python
Exploitation graphique
Exploitation graphique:
PCSI-MPSI – 9
Champs de vecteurs
Solution
Implémentation de la méthode d’Euler en Python
1.
On charge les modules numpy (ou scipy )pour les calculs numériques et
2.
3.
matplotlib.pyplot (ou pylab )pour les graphiques
La fonction (y; t) 7→ f (y; t) est définie au préalable.
On rappelle que y0 = y(t0 ) ∈ R n éventuellement
Version pédagogique
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
import numpy as np
import matplotlib . pyplot as mp
##
def Euler (f ,t0 , y0 , tf , N ):
h =( tf - t0 )/ N
# definition du pas fixe
# sinon h= t [i +1] - t [ i]
Ly , Lt =[ y0 ] ,[ t0 ]
# stockage du point initial
f o r k in range ( N ):
# A la limite Lt peut etre donne ailleurs
# prevoir que y0 peut etre vectoriel
y0 = y0 + h * np . array ( f (t0 , y0 ))
t0 += h
Ly . append ( y0 )
Lt . append ( t0 )
r e t u r n np . array ( Lt ) , np . array ( Ly )
##
A. Hassan@Champollion
Version optimisée
1
2
3
4
5
6
7
8
9
import numpy as np
def euler ( dXdt ,X0 , t ):
n = len ( t)
X = np . zeros (( len ( t ) ,)+ np . shape ( X0 ))
X [0]= X0
f o r i in range (n -1):
h =( t [ i +1] - t[ i ])
X [ i +1]= X [ i ]+ dXdt ( X [i ] , t[ i ])* h
return X
PCSI-MPSI – 10
Runge Kutta d’ordre 4 (RK4)
Soit (t; y(t)) un point courant, y(t + h) est une estimation de la
solution en t + h
Pour n = 4 ⇒
Méthode ou Runge Kutta d’ordre 4.
K1 = h. f (y(t); t);
1
2
1
2
K2 = h. f ( K1 + y(t); t + h)
1
1
K3 = h. f (y(t) + K2 ; t + h);
2
2
K4 = h f ( y ( t ) + K3 : t + h )
1
6
y(t + h) ≈ y(t) + (K1 + 2K2 + 2K3 + K4 )
Résolution des équation
différentielles ordinaires
(EDO)
Le Problème
Le Problème (suite)
Problème de Cauchy
Méthodes de résolution:
Ordres 1 et 2
Méthode d’Euler ou RK1
(Runge-Kutta d’ordre 1)
Implémentation de
l’algorithme d’Euler
Implémentation de la
méthode d’Euler en
Python
Runge Kutta d’ordre 4
(RK4)
Méthode de Runge-Kutta
d’ordre 4 en Python
Méthode de Runge-Kutta
d’ordre 4 en Python
Utilisation de la
commande odeint du
module
scipy.integrate
Il existe des méthodes à pas adaptatif (h n’est pas fixe).
A. Hassan@Champollion
Commande
odeint (suite)
Exercice: L’équation de
Van Der Pol (1924)
Implémentation en Python
Exploitation graphique
Exploitation graphique:
PCSI-MPSI – 11
Champs de vecteurs
Solution
Méthode de Runge-Kutta d’ordre 4 en Python
On charge au préalable le module numpy (ou scipy ). La fonction f est pointée
par dXdt .
1
import numpy as np
Voici un code de RK4 :
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
def rk4 ( dXdt ,X0 , t ): # Runge Kutta d ’ ordre 4
"""
Entrees :
t
: contient les points de la subdivision
X0
: les conditions initiales
dXdt
: la fonction connue
Sorties :
X
: Solutions approchees aux points de t """
# import numpy as np
n = len (t )
X = np . zeros (( len ( t ) ,)+ np . shape ( X0 )) # on construit une matrice de 0
X [0]= X0
# stocker les valeurs initiales
f o r i in range (n -1):
k1 =................
# k1 =f (y ,t )
h =................
k2 =................
k3 =................
k4 =................
X [i +1]= X [i ]+................
# on a deja divise h par 2
return (X)
A. Hassan@Champollion
PCSI-MPSI – 12
Méthode de Runge-Kutta d’ordre 4 en Python
On charge au préalable le module numpy (ou scipy )
1
import numpy as np
Voici un code de RK4 :
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
def rk4 ( dXdt ,X0 , t ): # Runge Kutta d ’ ordre 4
"""
Entrees :
t
: contient les points de la subdivision
X0
: les conditions initiales
dXdt
: la fonction connue
Sorties :
X
: Solutions approchees aux points de t """
# import numpy as np
n = len (t )
X = np . zeros (( len ( t ) ,)+ np . shape ( X0 )) # on construit une matrice de 0
X [0]= X0
# stocker les valeurs initiales
f o r i in range (n -1):
k1 = dXdt ( X[ i ] ,t [ i ])
# k1 = f(y , t)
h =( t [i +1] - t [ i ])/2
k2 = dX_dt (X [ i ]+ h * k1 , t [i ]+ h)
k3 = dX_dt (X [ i ]+ h * k2 , t [i ]+ h)
k4 = dX_dt (X [ i ]+2* h* k3 ,t [ i ]+2* h )
X [i +1]= X [i ]+ h /3*( k1 +2* k2 +2* k3 + k4 ) # on a deja divise h par 2
return (X)
A. Hassan@Champollion
PCSI-MPSI – 13
Utilisation de la commande odeint du module
scipy.integrate
Le module scipy.integrate fournit la commande odeint qui est un "mix" de
Runge Kutta d’ordre 4 avec des méthodes à pas adaptatif.
Syntaxe et utilisation :
1.
2.
3.
chargement du module adéquat :
from scipy.integrate import odeint
forme simplifiée :
odeint(func,CondInit,t)
forme complète :
odeint(func,
y0,
t,
args=(),
Dfun=None,
col_deriv=0, full_output=0, ml=None, mu=None,
rtol=None,
atol=None,
tcrit=None,
h0=0.0,
hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0,
mxordn=12, mxords=5, printmessg=0)
A. Hassan@Champollion
PCSI-MPSI – 14
Commande odeint (suite)
Forme simplifiée :
Ys=odeint(f, CondInit, t)
1.
2.
f pointe sur la fonction CONNUE f dans y′ (t) = f (y(t), t).
Attention : le premier argument de f est la fonction
vectorielle inconnue y.
CondInit pointe sur le vecteur des conditions initiales


y1 ( t0 )


y(t0 ) =  ... 
y n ( t0 )
3.
t pointe sur les points définissant le découpage de
l’intervalle de résolution I = [t0 ; t0 + T ]
La solution Ys contient les valeurs des solutions approchées
sur les différents points indiqués de I
A. Hassan@Champollion
Résolution des équation
différentielles ordinaires
(EDO)
Le Problème
Le Problème (suite)
Problème de Cauchy
Méthodes de résolution:
Ordres 1 et 2
Méthode d’Euler ou RK1
(Runge-Kutta d’ordre 1)
Implémentation de
l’algorithme d’Euler
Implémentation de la
méthode d’Euler en
Python
Runge Kutta d’ordre 4
(RK4)
Méthode de Runge-Kutta
d’ordre 4 en Python
Méthode de Runge-Kutta
d’ordre 4 en Python
Utilisation de la
commande odeint du
module
scipy.integrate
Commande
odeint (suite)
Exercice: L’équation de
Van Der Pol (1924)
Implémentation en Python
Exploitation graphique
Exploitation graphique:
PCSI-MPSI – 15
Champs de vecteurs
Solution
Exercice : L’équation de Van Der Pol (1924)
Il s’agit de l’équation différentielle
y′′ + f (y2 − 1)y′ + y = 0
( f ≥ 0)
Où f .(y2 − 1) représente un facteur d’amortissement non
linéaire ( f ∈ R + ).
1.
2.
3.
4.
Préciser à quoi correspond le cas f = 0 ? f 6= 0 et a 6= 0 ?.
Écrire l’équation de Van der Pol sous forme d’une équation
vectorielle de premier ordre.
Montrer que 0 est le seul point d’équilibre et qu’il n’est pas
stable.
Expliquer qualitativement le comportement lorsque y ≈ 0
ou bien lorsque y croît.
Pour les simulations on prend
kπ
(y0 ; y0′ ) ∈ {r (cos( kπ
)
;
sin
(
5
5 )) : k ∈ [[ 0, 9 ]] et r ∈ {1; 3}}
A. Hassan@Champollion
Résolution des équation
différentielles ordinaires
(EDO)
Le Problème
Le Problème (suite)
Problème de Cauchy
Méthodes de résolution:
Ordres 1 et 2
Méthode d’Euler ou RK1
(Runge-Kutta d’ordre 1)
Implémentation de
l’algorithme d’Euler
Implémentation de la
méthode d’Euler en
Python
Runge Kutta d’ordre 4
(RK4)
Méthode de Runge-Kutta
d’ordre 4 en Python
Méthode de Runge-Kutta
d’ordre 4 en Python
Utilisation de la
commande odeint du
module
scipy.integrate
Commande
odeint (suite)
Exercice: L’équation de
Van Der Pol (1924)
Implémentation en Python
Exploitation graphique
Exploitation graphique:
PCSI-MPSI – 16
Champs de vecteurs
Solution
Implémentation en Python
Préambule
1
2
3
4
# # modules a charger
from scipy import *
from matplotlib . pyplot import *
from scipy . integrate import odeint
Programmation
Van Der Pol :
1
2
3
4
5
6
7
def VdPol (......):
...........
...........
...........
...........
...........
#
Simulation
1
2
3
du système correspondant à l’équation de
#y′′ + f (y2 − 1)y′ + y = 0
et utilisation de odeint :
t = linspace (...........)
ci =[.... ,....]
Ys = odeint (..... , ..... , ..........)
A. Hassan@Champollion
Résolution des équation
différentielles ordinaires
(EDO)
Le Problème
Le Problème (suite)
Problème de Cauchy
Méthodes de résolution:
Ordres 1 et 2
Méthode d’Euler ou RK1
(Runge-Kutta d’ordre 1)
Implémentation de
l’algorithme d’Euler
Implémentation de la
méthode d’Euler en
Python
Runge Kutta d’ordre 4
(RK4)
Méthode de Runge-Kutta
d’ordre 4 en Python
Méthode de Runge-Kutta
d’ordre 4 en Python
Utilisation de la
commande odeint du
module
scipy.integrate
Commande
odeint (suite)
Exercice: L’équation de
Van Der Pol (1924)
Implémentation en Python
Exploitation graphique
Exploitation graphique:
PCSI-MPSI – 17
Champs de vecteurs
Solution
Exploitation graphique
y1 ( t )
Réponses Temporelles Si Y (t) =
solution de
y2 ( t )
Y ′ = f (Y, t). Les courbes temporelles correspondent aux
courbes (t; y1 (t)) et (t; y2 (t))
1
2
3
4
5
6
clf ()
# Tracer la reponse temporelle
plot (
,
, ’r ’) # premiere composante
plot (
,
, ’g ’) # seconde composante
grid ();
# grille
show ()
Portraits de phase : Le portrait de phase correspond la
trajectoire M (t) = (y1 (t); y2 (t))
1
2
CIs =[[ , ] ,[ , ] ,[ , ] ,[ , ]]
# ensemble de conditions initiales
f o r n in ................:
3
4
Ysol =...................
5
6
7
plot (................. .. ... ... ... )
show ()
A. Hassan@Champollion
Résolution des équation
différentielles ordinaires
(EDO)
Le Problème
Le Problème (suite)
Problème de Cauchy
Méthodes de résolution:
Ordres 1 et 2
Méthode d’Euler ou RK1
(Runge-Kutta d’ordre 1)
Implémentation de
l’algorithme d’Euler
Implémentation de la
méthode d’Euler en
Python
Runge Kutta d’ordre 4
(RK4)
Méthode de Runge-Kutta
d’ordre 4 en Python
Méthode de Runge-Kutta
d’ordre 4 en Python
Utilisation de la
commande odeint du
module
scipy.integrate
Commande
odeint (suite)
Exercice: L’équation de
Van Der Pol (1924)
Implémentation en Python
Exploitation graphique
Exploitation graphique:
PCSI-MPSI – 18
Champs de vecteurs
Solution
Exploitation graphique : Champs de vecteurs
On voudrait en chaque point M du plan faire apparaˆtre un
vecteur ~v tangent à la courbe intégrale passant par ce point (i.e.
~v colinéaire à (1; f (y, t)) dans notre cas)
1
2
3
4
5
6
7
8
9
10
11
12
# # Creation des champs de vecteurs
x = linspace ( -4 ,4 ,20)
# subdiviser l ’ intervalle de x [−4; 4]
y = linspace ( -4 ,4 ,20)
# subdiviser l ’ intervalle de y [−4; 4]
X1 , Y1 = meshgrid (x ,y )
# grille de noeuds en x , y
dX , dY = VdPol ([ X1 , Y1 ] ,0) # generer les vecteurs tangents
M = hypot ( dX , dY )
# normalisation
M [M ==0]=1.
# Normes d ’ eventuels vecteurs nuls
dX /= M
# remplaces par 1 avant division
dY /= M
quiver ( Y1 ,X1 ,dY , dX , M )
# generation du champs de vecteurs
show ()
##
Résolution des équation
différentielles ordinaires
(EDO)
Le Problème
Le Problème (suite)
Problème de Cauchy
Méthodes de résolution:
Ordres 1 et 2
Méthode d’Euler ou RK1
(Runge-Kutta d’ordre 1)
Implémentation de
l’algorithme d’Euler
Implémentation de la
méthode d’Euler en
Python
Runge Kutta d’ordre 4
(RK4)
Méthode de Runge-Kutta
d’ordre 4 en Python
Méthode de Runge-Kutta
d’ordre 4 en Python
Utilisation de la
commande odeint du
module
scipy.integrate
Commande
odeint (suite)
A. Hassan@Champollion
Exercice: L’équation de
Van Der Pol (1924)
Implémentation en Python
Exploitation graphique
Exploitation graphique:
PCSI-MPSI – 19
Champs de vecteurs
Solution
Solution
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
# # Resolution de l ’ equation de Van der Pol y′′ + f ∗ (y2 − 1)y′ + y = 0
def VdPol (y , t ): #
""" Attention d ’ abord y ( t) puis t """
f =1
v =y [0];
# y
vp = y [1];
# y′
vpp = f *(1 - v **2)* vp - v # y′′ = −y + f (1 − y2 )y′
r e t u r n array ([ vp , vpp ])
##
t = linspace (0 ,20 ,1000)
# on simule sur 1000 points de T = [0; 20]
# genere des conditions initiales
ci = array ([[0.5* cos ( k * pi /5) ,0.5* sin ( k* pi /5)] f o r k in range (10)])
ci2 = array ([[2.5* cos ( k* pi /5) ,2.5* sin (k * pi /5)] f o r k in range (10)])
#
f o r i in range (10):
Ys = odeint ( VdPol , ci [ i ] ,t )
Yys = odeint ( VdPol , ci2 [i ] , t)
Vp = Ys [: ,1]; V= Ys [: ,0]
# recuperer les reponses temporelles
plot ( Vp ,V , ’r ’)
# le portrait de phase ( rouge )
plot ( Yys [: ,1] , Yys [: ,0] , ’g ’) # encore des portraits de phases
plot ( ci [: ,0] , ci [: ,1] , ’o ’)
plot ( ci2 [: ,0] , ci2 [: ,1] , ’ob ’)
title (" Trajectoires de Van der Pol : y′′ + f (y2 − 1)y′ + y = 0" )
xlabel ( " Apparition du cycle limite " ); grid ()
savefig ( " G :\\ Azzam \\ MyPython \\ TpChampollion \\ VanDerPol . eps ")
show ()
A. Hassan@Champollion
PCSI-MPSI – 20
Portrait de phase
3
Trajectoires de Van der Pol:
Résolution des équation
différentielles ordinaires
(EDO)
Le Problème
Le Problème (suite)
Problème de Cauchy
Méthodes de résolution:
Ordres 1 et 2
Méthode d’Euler ou RK1
(Runge-Kutta d’ordre 1)
y′′+ f(y2 −1)y′ + y =0
2
1
Implémentation de
l’algorithme d’Euler
Implémentation de la
méthode d’Euler en
Python
Runge Kutta d’ordre 4
(RK4)
Méthode de Runge-Kutta
d’ordre 4 en Python
Méthode de Runge-Kutta
d’ordre 4 en Python
Utilisation de la
commande odeint du
0
−1
−2
−3
−4
A. Hassan@Champollion
module
scipy.integrate
−3
−2
−1
1
0
Apparition du cycle limite
2
3
4
Commande
odeint (suite)
Exercice: L’équation de
Van Der Pol (1924)
Implémentation en Python
Exploitation graphique
Exploitation graphique:
PCSI-MPSI – 21
Champs de vecteurs
Solution
Solution(suite)
Réponses temporelles :
1
2
3
4
5
6
7
8
9
10
11
subplot (221)
plot (t , Vp ); grid ();
title (" évolution de y′ (t) " )
subplot (222)
plot (t , V ); grid ();
title (" évolution de y(t) ")
subplot (212)
plot ( Vp , V ); grid ();
title (" portrait de phase correspondant ")
show ()
savefig ( " C :// Users // David // Desktop // Info // Python // Cours // analysenum // RepTempo . eps ")
A. Hassan@Champollion
PCSI-MPSI – 22
Solution : Graphiques(suite)
évolution de y (t)
′
3
3
2
2
1
1
0
0
−1
−1
−2
−2
−3
0
5
10
15
20
−3
0
évolution de y(t)
5
15
10
portrait de phase correspondant
3
Résolution des équation
différentielles ordinaires
(EDO)
Le Problème
Le Problème (suite)
Problème de Cauchy
Méthodes de résolution:
Ordres 1 et 2
Méthode d’Euler ou RK1
(Runge-Kutta d’ordre 1)
20
2
1
0
−1
module
scipy.integrate
−2
−3
−3
A. Hassan@Champollion
Implémentation de
l’algorithme d’Euler
Implémentation de la
méthode d’Euler en
Python
Runge Kutta d’ordre 4
(RK4)
Méthode de Runge-Kutta
d’ordre 4 en Python
Méthode de Runge-Kutta
d’ordre 4 en Python
Utilisation de la
commande odeint du
−2
−1
0
1
2
3
Commande
odeint (suite)
Exercice: L’équation de
Van Der Pol (1924)
Implémentation en Python
Exploitation graphique
Exploitation graphique:
PCSI-MPSI – 23
Champs de vecteurs
Solution
Solution : Champs de vecteurs
1
2
3
4
5
6
7
8
9
10
11
12
13
clf ()
x = linspace ( -4 ,4 ,30)
y = linspace ( -4 ,4 ,30)
X1 , Y1 = meshgrid (x ,y )
dX , dY = VdPol ([ X1 , Y1 ] ,0)
M = hypot ( dX , dY )
M [M ==0]=1.0
dX /= M
dY /= M
#
quiver ( Y1 ,X1 ,dY , dX , M )
show ()
savefig ( " C :// Users // David // Desktop // Info // Python // Cours // analysenum // ChampVdPol . eps " )
A. Hassan@Champollion
PCSI-MPSI – 24
4
3
2
1
0
−1
−2
−3
−4
−4
A. Hassan@Champollion
−3
−2
−1
0
1
2
3
4
PCSI-MPSI – 25
Solution : Superposition Portrait de phase/ champs de
vecteurs
4
3
2
1
0
−1
−2
−3
−4
−4
A. Hassan@Champollion
−3
−2
−1
0
1
2
3
4
PCSI-MPSI – 26
Exemple : Le pendule simple (non amorti)
Soit une masse m suspendu à un point fixe par un fil non extensible de longueur ℓ.
Les lois de la mécanique impliquent :
g
l
mℓ2 θ ′′ = −mgℓ sin(θ ) ⇐⇒ θ ′′ = − sin(θ )
1.
θ angle de l’écart avec la verticale en radians, θ ′ (t) (resp. θ ′′ (t)) vitesse (resp. accélération)angulaire
2. ℓ longueur en mètre.
3. g accélération de la pesanteur en m.s−2
4. À l’instant t = 0, θ = θ0 écart par rapport à la
verticale et θ0′ = 0 (vitesse initiale nulle)
bc
θ
bc
m~g
Question : Trouver la loi du mouvement t 7→ θ (t) ?
Pas de solution exprimable à l’aide des fonction usuelles (composées de exp() et
de ln() . . . ), (on démontre même que cela est impossible).
Remarque : Résolution pour des "petits" θ :
q
sin( θ )≈ θ
=======
⇒ θ (t) = θ0 cos(ωt) où ω = gℓ
A. Hassan@Champollion
PCSI-MPSI – 27
Traitement
On cherche une fonction y : R + −→ I × R où I = [−π; π ] où
θ position angulaire
θ
où
y=
θ ′ vitesse angulaire
θ′
#
"
′ ′
y2
θ
θ
′
g
g
= − sin(θ ) =
y =
− sin(y1 )
θ ′′
ℓ
ℓ
Donc
f : I × R −→ R × R
g
( x1 ; x2 ) 7→ ( x2 ; − sin( x1 ))
ℓ
L’équation devient (Problème de Cauchy) :
y′ = f (y)
avec y(0) =
θ0
0
Équation différentielle autonome(stationnaire) : f ne dépend
pas de t explicitement.
A. Hassan@Champollion
Résolution des équation
différentielles ordinaires
(EDO)
Le Problème
Le Problème (suite)
Problème de Cauchy
Méthodes de résolution:
Ordres 1 et 2
Méthode d’Euler ou RK1
(Runge-Kutta d’ordre 1)
Implémentation de
l’algorithme d’Euler
Implémentation de la
méthode d’Euler en
Python
Runge Kutta d’ordre 4
(RK4)
Méthode de Runge-Kutta
d’ordre 4 en Python
Méthode de Runge-Kutta
d’ordre 4 en Python
Utilisation de la
commande odeint du
module
scipy.integrate
Commande
odeint (suite)
Exercice: L’équation de
Van Der Pol (1924)
Implémentation en Python
Exploitation graphique
Exploitation graphique:
PCSI-MPSI – 28
Champs de vecteurs
Solution
Exemples : Résolution de y′′ + y = 0 (petite oscillation du
pendule)
On dispose d’une solution exacte de y′′ + y = 0 (n’est-ce pas ?) où
(y(0); y′ (0)) = ( a, b) :
y(t) = a cos(t) + b sin(t)
1
2
3
4
5
6
7
# # Un exemple : y ’ ’+y =0
def f (y , t ): # definition de la fonction vectorielle
""" y (1)= ypp ,
y (2)= yp """
r e t u r n ( np . array ([ - y [1] , y [0] ]))
#
t = np . linspace (0 ,30 ,1000) # découpage de l ’ intervalle [0 ,10] en 1000
8
9
10
11
12
13
14
15
16
17
18
19
20
# utilisation des méthodes d ’ Euler et RK4 programmées
# cond . init . (y′ (0), y(0)) = (0; 1)
ci = np . array ([0 ,1])
# solution par Euler
sol1 = euler (f ,ci , t)
# solution par RK4
sol2 = rk4 (f , ci , t )
# solution exacte de y ’ ’+ y =0
solexact = np . cos ( t) # y(t) = cos(t)
# utilisation de la methode odeint de numpy . integrate
sol3 = sint . odeint (f , ci , t )
#
A. Hassan@Champollion
PCSI-MPSI – 29
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# Exploitation graphique
py . subplot (211)
py . title ( " Comparaison : Euler ( rouge ) , solution exacte ( noir ) de y′′ + y = 0" )
py . plot (t , sol1 [: ,1] , color = ’ red ’)
py . plot (t , solexact , color = ’ black ’)
#
py . subplot (212)
py . title ( " Comparaison : RK4 ( rouge ) , solution exacte de y′′ + y = 0" )
py . plot (t , sol2 [: ,1] , color = ’ red ’)
py . plot (t , solexact , color = ’ black ’)
py . grid ( True )
py . savefig ( " edo1 . eps ") # sauvegarder l ’ image dans un fichier , SI JE VEUX
py . subplot (111)
py . title ( " Comparaison : numpy . odeint ( rouge ) , solution exacte de y′′ + y = 0")
py . plot (t , sol3 [: ,1] , color = ’ red ’)
py . plot (t , solexact , color = ’ black ’)
py . grid ( True )
py . savefig ( " edo2 . eps ")
A. Hassan@Champollion
PCSI-MPSI – 30
Résultats graphiques
Résolution des équation
différentielles ordinaires
(EDO)
Le Problème
Le Problème (suite)
Problème de Cauchy
Méthodes de résolution:
Ordres 1 et 2
Méthode d’Euler ou RK1
(Runge-Kutta d’ordre 1)
Comparaison: Euler(rouge), solution exacte(noir) de y + y =0
1.5
′′
1.0
0.5
0.0
−0.5
−1.0
−1.5
−2.0
0
1.0
5
10
15
20
25
Comparaison: RK4(rouge), solution exacte de y + y =0
′′
30
0.5
0.0
module
scipy.integrate
−0.5
−1.0
0
A. Hassan@Champollion
Implémentation de
l’algorithme d’Euler
Implémentation de la
méthode d’Euler en
Python
Runge Kutta d’ordre 4
(RK4)
Méthode de Runge-Kutta
d’ordre 4 en Python
Méthode de Runge-Kutta
d’ordre 4 en Python
Utilisation de la
commande odeint du
5
10
15
20
25
30
Commande
odeint (suite)
Exercice: L’équation de
Van Der Pol (1924)
Implémentation en Python
Exploitation graphique
Exploitation graphique:
PCSI-MPSI – 31
Champs de vecteurs
Solution
Résultats graphiques
Comparaison: numpy.odeint(rouge), solution exacte de y + y =0
1.0
′′
0.5
Résolution des équation
différentielles ordinaires
(EDO)
Le Problème
Le Problème (suite)
Problème de Cauchy
Méthodes de résolution:
Ordres 1 et 2
Méthode d’Euler ou RK1
(Runge-Kutta d’ordre 1)
Implémentation de
l’algorithme d’Euler
Implémentation de la
méthode d’Euler en
Python
Runge Kutta d’ordre 4
(RK4)
Méthode de Runge-Kutta
d’ordre 4 en Python
Méthode de Runge-Kutta
d’ordre 4 en Python
Utilisation de la
commande odeint du
0.0
−0.5
module
scipy.integrate
−1.0
0
A. Hassan@Champollion
5
10
15
20
25
30
Commande
odeint (suite)
Exercice: L’équation de
Van Der Pol (1924)
Implémentation en Python
Exploitation graphique
Exploitation graphique:
PCSI-MPSI – 32
Champs de vecteurs
Solution
Pendule simple : y′′ + ω2 sin(y) = 0
Cette fois on ne fait pas l’approximation sin(y) ≈ y.
Pas de solution exacte simple, la comparaison avec la solution "petits angles"
donne :
1
2
3
# # Resolution de y ’ ’+ sin ( y )=0
def f (y , t ):
r e t u r n ( np . array ([ - np . sin (y [1]) , y [0] ]))
1.5
Euler(rouge)de
y
y) = 0,
′′
+ sin(
solution petits angles:y
′′
+
y =0
1.0
0.5
0.0
−0.5
−1.0
−1.5
0
1.0
5
RK4(rouge) de
y
′′
10
y) = 0,
+ 2sin(
15
20
25
solution petits angles:y
′′
+
y =0
30
0.5
0.0
−0.5
−1.0
0
A. Hassan@Champollion
5
10
15
20
25
30
PCSI-MPSI – 33
Résultats graphiques
1.0
odeint(rouge) de
y
′′
y) = 0,
+ sin(
solution petits angles:y
′′
+
y =0
0.5
0.0
−0.5
−1.0
0
A. Hassan@Champollion
5
10
15
20
25
30
PCSI-MPSI – 34
Construction du portrait de phases
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
t = np . linspace (0 ,7 ,300)
# On stocke pour plusieurs cond init .
courbes =[\
([0 ,1] , t ) ,\
([0 ,2] , t ) ,\
([0 ,4] , t ) ,\
([0 ,5] , t )]
#
f o r i in range ( len ( courbes )):
sol = sint . odeint (f , courbes [ i ][0] , courbes [i ][1])
py . plot ( sol [: ,0] , sol [: ,1]) # portrait de phase
# dessin du champs de vecteurs
# creation d ’ une grille
x = np . linspace ( -5 ,5 ,20)
y = np . linspace ( -5 ,5 ,20)
X ,Y = py . meshgrid (x , y )
# calcul des vecteurs du champ
DX , DY =f ([ X , Y ] ,0)
# normalisation
N = py . hypot ( DX , DY )
DX /= N ; DY /= N
# dessin du champ le dernier paramètre permet \
# de coloriser en fonction de la norme
py . quiver (X ,Y , DX ,DY , N)
py . title ( " champs et portrait de phase correspondant a y′′ + y = 0" )
py . grid ()
py . savefig ( " phase1 . eps " )
A. Hassan@Champollion
PCSI-MPSI – 35
Portrait de phase de y′′ + y = 0 : Visualisation
6
champs et portrait de phase correspondant a
y + y =0
′′
4
2
0
−2
−4
−6
−6
A. Hassan@Champollion
−4
−2
0
2
4
6
PCSI-MPSI – 36
Portrait de phase de y′′ + sin(y) = 0 : Visualisation
champ correspondant à l'équation du pendule:
y
′′
y) = 0
+ sin(
4
2
0
−2
−4
−4
A. Hassan@Champollion
−2
0
2
4
PCSI-MPSI – 37
Attracteur de Lorenz(vers 1963)
Les équations de Lorentz, basées sur la mécanique des fluides (équations de
Navier-Stokes), modélisent entre autre, l’évolution atmosphérique et expliquent
(en partie= les difficultés à rencontrées lors des prévisions météorologiques :
x ′ (t) = σ(y(t) − x (t))
y′ (t) = ̺x (t) − y(t) − x (t)z(t)
z′ (t) = x (t)y(t) − βz(t)
σ ≈ 10, β ≈ 2.66 ̺ ≈ 30
Où σ
x (t) y(t)
z(t)
A. Hassan@Champollion
: (dépend de la viscosité et la conductivité thermique)
: est proportionnel au taux de convection.
: des gradients de températures horizontal et vertical
PCSI-MPSI – 38
Résolutions des équations de Lorenz avec Python
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
import matplotlib . pyplot as mp
import numpy as np
from scipy . integrate import odeint
from mpl_toolkits . mplot3d import Axes3D # Module pour graphique 3 D
# # Attracteur de Lorentz
sigma , rho , b = 10 ,30 ,2 # valeurs des parametres
fig = mp . figure ()
ax = fig . gca ( projection = ’3d ’)
# v =[ x ,y , z] dans le modele de Lorenz
def fLorenz ( v , t ):
r e t u r n ([ sigma *( v [1] - v [0]) , rho * v [0] - v [1] - v [0]* v [2] ,\
v [0]* v [1] - b * v [2]])
#
def OrbiteLorenz ( ci ,n , T ):
t = np . linspace (0 ,T , n +1)
f o r y0 in ci :
values = odeint ( fLorenz , y0 , t )
ax . plot ( values [: ,0] , values [: ,1] , values [: ,2] , label = str ( y0 ))
OrbiteLorenz ([[20 ,0 ,0] ,[10 ,0 ,0] ,[0.1 ,0 ,0]] ,1000 ,10)
ax . legend ()
mp . savefig ( " Lorenz . eps " )
mp . show ()
A. Hassan@Champollion
Résolution des équation
différentielles ordinaires
(EDO)
Le Problème
Le Problème (suite)
Problème de Cauchy
Méthodes de résolution:
Ordres 1 et 2
Méthode d’Euler ou RK1
(Runge-Kutta d’ordre 1)
Implémentation de
l’algorithme d’Euler
Implémentation de la
méthode d’Euler en
Python
Runge Kutta d’ordre 4
(RK4)
Méthode de Runge-Kutta
d’ordre 4 en Python
Méthode de Runge-Kutta
d’ordre 4 en Python
Utilisation de la
commande odeint du
module
scipy.integrate
Commande
odeint (suite)
Exercice: L’équation de
Van Der Pol (1924)
Implémentation en Python
Exploitation graphique
Exploitation graphique:
PCSI-MPSI – 39
Champs de vecteurs
Solution
Attracteur de Lorenz
[20, 0, 0]
[10, 0, 0]
[0.1, 0, 0]
60
50
40
30
20
10
−20−15
−10 −5
0
−10
5 10
−20
15 20
−30
25
0
10
20
0
30
Impossible de prévoir à long terme l’état du système.
A. Hassan@Champollion
PCSI-MPSI – 40
Le module NUMPY : quelques commandes
Vecteurs
1
2
3
4
5
vl = array ([ a ,b ,..])
vc = arry ([[ a ] ,[ b ] ,..])
vl [ i ]; vc [ j ]
vl . size
S = arange (d ,f , i )
#
#
#
#
#
vecteur ligne
vecteur colonne
elemt : vli et vc j
longueur de vl
S =[d , d +i , d +2i ,...] ( f 6∈ S)
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
remplissage par lignes
mij
renvoie le format de mat
extrait la ligne i de mat
extrait la colonne j de mat
0n × p
mij = 1 si i = j, 0 sinon
matrice des uns ( matrice d ’ Attila )
no comment
re - dimensionner : mat ∈ Mm×q (K )
transformer matrice en tableau
transformer tableau -> matrice
det(mat) avec le module numpy.linalg
mat−1 avec le module numpy.linalg
valeurs et vecteurs propres de mat avec le module num
vecteurs propres avec le module numpy.linalg
Matrices
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
mat =[[ a1 , a2 ,..] ,[ b1 , b2 ,..]]
mat [i , j ]
mat . shape
mat [i ,:]
mat [: , j ]
zeros ([ n , p ]); zeros (n , p)
m = eye (n ,p )
ones ([ n ,p ])
transpose ( mat )
reshape ( mat ,( m ,q ))
asarray ( mat )
asmatrix ( m )
det ( mat )
inv ( mat )
eigvals ( mat )
eig ( mat )
A. Hassan@Champollion
PCSI-MPSI – 41
Commandes NUMPY (suite)
Polynômes
1
2
3
4
5
6
P = poly1d ([ a ,b ,..])
Pf = poly1d ([ a ,b ,..] , True )
P . order
P . roots ; roots (P )
P . coeffs
P ( x ); polyval (P , x )
A. Hassan@Champollion
# P ( x) = axn + bxn−1 + · · ·
# Pf ( x) = ( x − a)( x − b) · · ·
# renvoie les racines complexes ( approchees )
# renvoie les coefficients de P
# renvoie P ( x) selon Horner
PCSI-MPSI – 42
Opération sur les matrices (resp. vecteurs)
1.
2.
3.
4.
5.
Combinaisons linéaires
Produit de matrices
inverse de matrices
Produit termes à termes
Produit tensoriel
A. Hassan@Champollion
PCSI-MPSI – 43