Sorbonne Université, Master 1 Mathématiques et applications 2024–2025

4MA006 & 4MA106 Fondements des méthodes numériques

Projet TP 1

Notation. On rappelle que votre note de TP correspond à 30% de votre note de 4MA006
(4MA106) et que cette note de TP est incompressible, c’est à dire sans rattrapage. Votre
projet est individuel et il sera évalué comme suit :
• rapport et code : clarté, lisibilité, organisation.
• soutenance orale : aisance à utiliser le code, pertinence des explications et réponses aux
questions, recul sur le lien entre résultats théoriques et observations numériques. La
soutenance orale représentera une très moitié de la note finale.
Vous êtes bien sûr autorisé à discuter du projet entre vous et avec l’équipe enseignante de l’UE
mais par individuel on entend que toutes similitudes trop importantes pour être considérées
comme le fruit d’un heureux hasard, et ceci que ce soit au niveau du code et/ou du rapport,
seront considérées comme de la triche et donc sanctionnées sévèrement au niveau de la note.

Le code. Le code doit être envoyé sous la forme d’un ou plusieurs fichiers et doit être
développé en Python (extention .py ou .ipynb ). Votre code doit être agréable à lire (noms de
variables explicites, indentation, un minimum de commentaires,...). Pour éviter tout problème
d’encodage, ne pas utiliser d’accent dans vos commentaires.

Le rapport. Le rapport doit être retourné en format .pdf. Dans le cas d’un rapport pdf, de
préférence utiliser LATEX mais vous pouvez également envoyer des notes manuscrites scannées
(soigneuses).

Deadline. Le nom de tous vos fichiers doit débuter par NOM prenom. Vos fichiers, code et
rapport, seront à déposer sur la page Moodle de 4M006 (4M106) avant le 30/12/2024.

Soutenance. La soutenance orale aura lieu la semaine des examens. Un planning sera
disponible sur Moodle en temps voulu. Cette soutenance consistera à présenter et commenter
votre travail. Vous disposerez de 10 minutes pour ce faire, puis s’en suivra une séance de
questions. Il vous est fortement conseillé de vous entraı̂ner à cet exercice avant votre passage.
Tout retard, ou non respect des consignes décrites ici, sera sanctionné dans votre
note de TP.

Le sujet comporte des parties indépendantes. L’important n’est pas de le finir, l’évaluation
portera plus sur la qualité du travail rendu que sa quantité. Pour toute question ou commen-
taire sur ce projet, il ne faut pas hésiter à prendre contact via Moodle ou par mail :

[email protected]

1
1 Approximation de l’équation de transport à vitesse con-
stante
Soit un temps final T > 0, nous cherchons une approximation de, ū : R × (0, T ) → R, solution
du problème suivant,
(
∂t ū(x, t) + a ∂x ū(x, t) = 0, ∀(x, t) ∈ R × (0, T ) ,
(P)
ū(x, 0) = uini (x), ∀x ∈ R .

On rappelle que la solution de (P) est donnée par ū(x, t) = uini (x − at), ∀(x, t) ∈ R × (0, T ).
Afin de travailler sur un domaine d’espace de dimension finie, on prendra une condition initiale
uini 1-périodique.

Exercice 1 Montrer que, puisque uini est 1-périodique, alors ū est également 1-périodique
en espace, i.e. ū(x + 1, t) = ū(x, t), ∀(x, t) ∈ R × (0, T ).

Exercice 2 On considère la donnée initiale, uini (x) = sin(2πx) , ∀x ∈ R. Vérifier qu’elle est
bien 1-périodique.

On définit une solution numérique discrète de l’équation comme une suite de valeurs unj à
double indice approchant la solution exacte ū(x, t) au sens où unj ≈ ū(xj , tn ) avec xj = j∆x
et tn = n∆t, pour j ∈ Z et n ∈ N. On nomme ∆x le pas d’espace et ∆t le pas de temps. On
∆t
pose λ = ∆x .

Commençons par quelques exemples de schémas explicites pour l’équation de transport.

1. Schéma centré. Ce schéma est obtenu en utilisant une approximation centrée de la
dérivée spatiale.
λa n
un+1 = unj − (u − unj−1 ).
j
2 j+1
2. Schéma décentré à gauche. Il est obtenu en utilisant une approximation décentrée
(à gauche) de la dérivée spatiale.
ujn+1 = unj − λa(unj − unj−1 ).

3. Schéma décentré à droite. Il est obtenu en utilisant une approximation décentrée

(à droite) de la dérivée spatiale.
un+1
j = unj − λa(unj+1 − unj ).

4. Schéma de Lax-Friedrichs. On reconnaı̂t une modification du schéma centré dans

lequel unj est remplacé par la moyenne de unj−1 et unj+1 .
1 λa n
un+1 = (unj+1 + unj−1 ) − (u − unj−1 ).
j
2 2 j+1
5. Schéma de Lax-Wendroff. On reconnaı̂t une modification du schéma centré dans
lequel on a rajouté un terme étrange qui contient une discrétisation de la dérivée seconde
en espace alors qu’il n’y a pas de dérivée seconde dans l’EDP considérée.
λa n λ 2 a2 n
un+1
j = unj − (uj+1 − unj−1 ) + (uj+1 − 2unj + unj−1 ).
2 2

2
On peut les écrire sous la forme
+L
X
un+1
j = H(unj−L , · · · , unj+L ) = cl unj+l , (1)
l=−L

où on a introduit une fonction H appelée “solution discrète”. Tous les schémas présentés sont
des schémas à 3 points, c’est-à-dire qu’ils n’utilisent qu’au plus, les valeurs aux trois points
xj−1 , xj et xj+1 pour avancer en temps (sous la forme , L = 1).

Exercice 3 On va donc implémenter le schéma de sorte à ce qu’il soit lui aussi 1-périodique.
La discrétisation en espace doit donc couvrir au minimum un intervalle de la forme [S, 1 + S[,
avec S ∈ R. Le plus naturel est de prendre S = 0 et donc l’intervalle [0, 1[.
1. On se fixe J ∈ N∗ et on pose ∆x = J1 . Définir une suite de points xj = j∆x adaptée à
la discrétisation de l’intervalle [0, 1[. On pose λ = ∆t/∆x.
(Conseil : Vous pouvez compléter une fonction espaceDiscretise dont les entrées et
sorties :
Input : J, λ.
Output : ∆x, ∆t, (xj )j , M .)

2. Coder la fonction qui s’appelle schema en utilisant la fonction discrète H dont voici
les entrées/sorties :
Input : ∆x, ∆t, (xj )j , M , c−1 , c0 , c1
Output : (uM M
j )j , (ūj )j
Dans le code, tM ≤ T correspond au temps final de simulation et ūM
j = ū(xj , tM ), ∀j.

3. Déterminer c−1 , c0 et c1 tel que chaque schéma se réécrive sous la forme,

un+1
j = c−1 unj−1 + c0 unj + c1 unj+1 ,

et tracer la solution analytique et numérique. Pour la discrétisation en espace déterminée

dans la première étape de cet exercice, réfléchir à la prise en compte de la périodicité en
espace au niveau du schéma à chaque itération en temps tn → tn+1 pour n ≥ 0. Dans
ce cas, on prend T = 0.75, a = 1, J = 20, λ = 0.8.

2 Analyse des schémas

L’analyse numérique des schémas a pour objectif d’étudier leur convergence, dans une norme
appropriée, et lorsqu’ils convergent, d’estimer l’erreur, calculée (souvent) dans la même norme,
entre la solution exacte et la solution approchée. Comme nous l’avons déjà vu pour l’équation
de la chaleur, deux propriétés importantes entrent en jeu : la stabilité et la consistance, dont
on déduit la convergence.

Commençons par l’étude de consistance. Comme pour l’équation de la chaleur, on adopte la

notation suivante.
Définition 2.1 On appelle erreur de consistance du schéma (1) au point xj et à l’instant tn ,
le réel
1
κnj = (u(xj , tn+1 ) − H (u(xj−L , tn ), . . . , u(xj+L , tn ))) .
∆t

3
L’erreur de consistance du schéma à l’instant tn est le vecteur Khn dont les composantes sont
les κnj , et qui vérifie donc
n+1 n
Uh = QU h + ∆tKhn . (2)

Définition 2.2 On dit qu’un schéma est stable pour une norme vectorielle k·k si et seulement
si il existe une constante C0 (qui peut dépendre de T ), telle que

sup kQn k ≤ C0 .
n∆t≤T

Exercice 4 Nous allons étudier la stabilité du schéma. Calculer numériquement la norme

matricielle pour chaque schéma. Dans ce cas, on prend T = 0.75, a = 1, J = 20, λ = 0.8.
Qu’observez-vous? Ensuite, Trouver la condition en écrivant à la main pour chaque schéma
(si possible) afin de déterminer sa stabilité.

Exercice 5 Nous devons étudier la convergence pour chaque schéma. On calcul la vitesse
à laquelle l’erreur décroit lors qu’on diminue ∆x (en aumentant J, le nombre des points).
On dit que la méthode est d’ordre p ∈ R+ s’il existe une constante C > 0, independente de
n n k = O((∆x)p ). On utilisera les normes Lq discrète. Pour différents
∆x, telle que kU ∆x − U∆x
n n k n n
valeurs de ∆x, calculer les erreurs kU ∆x − U∆x 2,M and kU ∆x − U∆x k∞,M . Representer les sur
une figure, en échelle log–log, en fonction de ∆x et déterminer les ordres. On liste ci-dessous
les expressions des normes Lq discrète
(∞)

∆t,∆x = max kU ∆x − U∆x k∞,M
n≥0
(2)

∆t,∆x = max kU ∆x − U∆x k2,M
n≥0

n = (un ) et Ū n = (ū(x , t )) . Dans ce cas, on prend T = 0.75, a = 1,

où U∆x j j ∆x j n j
J = [25, 50, 100, 200], λ = 0.8. Qu’en observez-vous?

(Conseil : Vous pouvez modifier la fonction schema dont les nouvelles entrées/sorties :
Input : ∆x, ∆t, (xj )j , M , c−1 , c0 , c1
(2) (∞)
Output : (uM M
j )j , (ūj )j , ∆t,∆x , ∆t,∆x .
Une fois vous avez la nouvelle fonction schema, vous pouvez obtenir les erreurs avec les
différents valeurs de ∆x. )