ENSAM-MEKNES 3 EME ANNEE

TP 2 METHODES NUMERIQUES
POUR INGENIEUR 2 2024-2025

Résolution numérique des EDOs

Exercice 1 Euler explicite

On cherche à mettre en œuvre la méthode d’Euler explicite pour la résolution numérique du
problème de Cauchy
y 0 (t) = f (t, y(t)), t ∈ [0, T ]
(

y(0) = y0 .
Elle consiste, pour un pas h donné, à construire une suite d’approximations (yn ) de la solution
y aux temps tn = nh, par la formule
yn+1 = yn + hf (tn , yn ).
Par la suite, on utilise l’exemple modèle y 0 = t − ty avec la condition initiale y(0) = 2. La
−t2
solution exacte de cet exemple est donnée par y(t) = 1 + e 2 .
A la fin de cette série de TP vous trouvez le programme principal, le programme du schéma
d’Euler explicite et le programme de la fonction f .
1. Tester le programme et interpréter la figure affichée
2. Programmer la solution exacte, puis calculer l’erreur d’approximation;
3. Tracer les deux solutions sur la même figure;
4. Effectuer des simulations pour différentes valeurs du pas h (h = T /N avec N = 10 : 10 :
1000);
5. Tracer la courbe du log des erreurs d’approximation en fonction du log des pas h;
6. Quel est l’ordre obtenue numériquement du schéma d’Euler explicite.
L’ordre numérique est la pente de la droite de la régression linéaire des points.

Exercice 2 Euler point milieu

On utilise ici l’exemple modèle précédent y 0 (t) = t − ty(t), y(0) = 2. La solution exacte est
−t2
donnée par y(t) = 1 + e 2 .
Le schéma d’Euler point milieu consiste à construire une suite d’approximations (yn ) de la
solution y de la façon suivante :
yn+1/2 = yn + h2 f (tn , yn )
yn+1 = yn + hf (tn + h2 , yn+1/2 ).
1. Ecrire une fonction "EulerPointMilieu" qui permet de calculer la solution approchée yn
par la méthode d’Euler point milieu.
2. Tracer les deux solutions approchées (Euler explicite et point milieu) et la solution exacte
sur la même figure
3. Déterminer numériquement l’ordre de la méthode d’Euler point milieu.
4. Programmer le schéma semi-implicite et le schéma de Heun.
5. Programmer le schéma d’Euler implicite.

1
 ENSAM-MEKNES 3 EME ANNEE
TP 2 METHODES NUMERIQUES
POUR INGENIEUR 2 2024-2025

Exercice 3 Runge Kutta : Application aux oscillations du pendule

On considère une pendule simple : la masse

ponctuelle m située au point M est soumise
à son poids P , et à la tension T du câble de
longueur L qui la relie au point fixe O. La
position de M est repérée par l’angle θ(t) en-
tre la verticale câble. En projection sur l’axe
u perpendiculaire à OM ~ , le théorème de la ré-
sultante dynamique donne :
00
mLθ (t) = −mg sin(θ(t)) (1)
00
où (Lθ ) est l’accélération linéaire de la masse

m dans la direction u, et − mg sin(θ) est
la force de rappel qui ramène m vers sa po-
sition d’équilibre. On en tire immédiatement
l’équation différentielle exacte:
00 g
θ =− sin(θ) (2)
L
1. Résolution de l’équation linéarisée
Pour les petits mouvements (θ ≈ 0), on linéarise l’équation différentielle en posant :
sin(θ) ≈ θ.
q 
(a) Écrivez l’équation linéarisée et vérifiez que θ(t) = A cos g
L
t+φ est une solution
dans ce cas.
• Pour g = 9.81m\s2 et L = 1m, quelles sont la pulsation, la fréquence et la
période du mouvement ? (Ce résultat classique correspond au fonctionnement
de l’horloge comtoise, dont le balancier de 1m de longueur “bat la seconde”).

(b) On remarque que, en remplaçant sin(θ) par θ, la linéarisation surestime la force de

rappel du pendule.
• Le mouvement calculé est-il plus ou moins rapide que le mouvement réel?
• La période calculée est-elle plus ou moins grande que la période réelle?

π 0
(c) On prend pour conditions initiales θ(0) = et θ (0) = 0. L’hypothèse des petits
2
mouvements n’est pas réalisée, et l’équation linéarisée donnera des résultats large-
ment entachés d’erreur.

• Déterminez les constantes d’intégration A et φ.

• Tracez θ(t), pour 0 6 t 6 4s, avec un pas de 0.04s. On peut suivre les étapes
suivantes :
i. Définissez une discrétisation T de l’intervalle [0, 4] avec un pas h = 0.04.
ii. Calculer le vecteur theta, composé des valeurs θ(t), ∀t ∈ T .
iii. Tracez theta en fonction de T.

2
 ENSAM-MEKNES 3 EME ANNEE
TP 2 METHODES NUMERIQUES
POUR INGENIEUR 2 2024-2025

2. Résolution numérique de l’équation exacte, par la méthode d’Euler

On rappelle le principe de la méthode d’Euler explicite :

yn+1 = yn + hF (tn , yn ). (3)

L’équation différentielle exacte du pendule peut se mettre sous la forme de 2 équations :

0
θ = ω


0 g (4)
 ω = − sin(θ).
L
Ceci peut s’exprimer sous la forme :
0
Y = F (t, Y ), (5)
!
θ
où Y est le vecteur défini par Y = , et F (t, Y ) est une fonction vectorielle.
ω

(a) Expliciter la fonction F (t, Y ).

(b) Programmer cette fonction avec Python sous le nom Fonct(t,Y)
(c) Programmez la méthode Euler explicite sous forme d’une fonction EulerExp(f,y0,a,b,h).
Cette fonction renvoie la discrétisation du temps T et le vecteur des solutions ap-
prochées Y .

(d) Résolvez ce problème :

!
θ(0)
i. Définissez le vecteur des conditions initiales Y 0 = et le pas de temps
ω(0)
h = 0.04s.
ii. Initialisez un tableau Ye à deux lignes et Npas=101 colonnes pour stocker les
valeurs approchées successives de Y (t).
iii. Calculer ces valeurs, de proche en proche, à l’aide de la fonction EulerExp
iv. Représentez les valeurs approchées successives de θ(t) sur la même figure que
les solutions du problème linéarisé.

(e) Que pensez-vous de l’amplitude et de la période du mouvement calculé par cette

méthode?
3. Résolution numérique de l’équation exacte, par la méthode de
Runge-Kutta
On rappelle le principe de la méthode de Runge-Kutta classique d’ordre 4 :
= F (tn , yn )


 k1
= F (tn + h/2, yn + 12 k1 h)

k2





= F (tn + h/2, yn + 21 k2 h)


k3
(6)
 k
 4

 = F (tn + h, yn + k3 h)
k1 + 2k2 + 2k3 + k4


yn+1 = yn + h



6


(a) Définir une fonction RungeKutta qui résout le problème de Cauchy par le schéma de
Runge-Kutta. La fonction a le même prototype que EulerExp.

3
 ENSAM-MEKNES 3 EME ANNEE
TP 2 METHODES NUMERIQUES
POUR INGENIEUR 2 2024-2025

(b) Calculer l’approximation de θ(t) par la méthode de Runge-Kutta.

(c) Représenter les valeurs approchées successives de θ(t) sur la même figure que les
tracés précédents.

4. Résolution numérique de l’équation exacte, avec le solveur odeint

On considère la fonction odeint de la librairie "scipy.integrate". Cette fonction per-
met de résoudre numériquement une EDO (vous pouvez taper help(odeint) pour plus
d’information sur cette fonction).

(a) Résolvez le problème en utilisant le solveur odeint, pour t ∈ [0, 4], avec les conditions
initiales définies précédemment.
(b) Tracez les valeurs de θ(t) calculées.
(c) Comparez avec les tracés précédents.

4
 ENSAM-MEKNES 3 EME ANNEE
TP 2 METHODES NUMERIQUES
POUR INGENIEUR 2 2024-2025

________________________________
# Programme principal
import numpy as np
import matplotlib.pyplot as plt
T=2.
t0=0
u0=2.
N=10
h=T/N
t,u=SchemaExp(fonct,u0,t0,T,N)
# ue = solext(t) # doit etre programme
# er=np.linalg.norm((ue-u),ord=np.inf)
plt.plot(t,u,color=’k’,ls=’--’,label="solution approchée")
# plt.plot(t,ue,color=’b’,ls=’:’,label="solution exacte")
plt.xlabel(’time’)
plt.ylabel(’solution u’)
plt.title(’Solution approchee par Euler explicite’)
plt.legend()
plt.show()
________________________________

________________________________
#Schema Euler Explicite : fonction "SchemaExp"
#
def SchemaExp(f,u0,t0,T,N):
h=T/N
t=np.linspace(t0,T,N+1)
u=np.zeros(N+1)
u[0]=u0
for i in range(N):
u[i+1]=u[i]+h*f(t[i],u[i])
return t,u
______________________

________________________________
#programme de la fonction f
#
def fonct(t,u):
return t-t*u
________________________________

________________________________
# Calcul de l’ordre numerique
# lh = liste des pas
# ler : liste des erreurs correspondantes à lh
print(np.polyfit(lh,ler,1))
________________________________