M1 CSA MSNU – Modélisation et simulation numérique 2014/15
Feuille de TP 3
Exercice Méthode de tir
On continue à étudier le problème aux limites
(1) −u00 (x) + c(x)u(x) = f (x) (x ∈ ]0, 1[ ) ; u(0) = u0 , u(1) = u1 .
Dans la méthode de tir, on résout à la place de (1) un problème de Cauchy
(2) −u00 (x) + c(x)u(x) = f (x) (x ∈ ]0, 1[ ) ; u(0) = u0 , u0 (0) = v ,
avec une inconnue supplémentaire v ∈ R . Le but est de choisir le paramètre v de façon à
ce que l’autre condition aux limites u(1) = u1 soit satisfaite. En notant F l’application
qui à v associe la valeur u(1) de la solution de (2), on cherche donc une solution v de
l’équation
(3) F (v) − u1 = 0.
Dans la suite, on écrira une fonction Matlab FF qui réalise (une approximation de) la
fonction F , et on résoudra l’équation (3) par la méthode de dichotomie (bissection).
1. Pour la solution du problème de Cauchy (2), on utilisera la fonction ode23 de Mat-
lab. Cette fonction calcule des solutions approchées d’un problème de Cauchy pour un
système d’ordre 1. On doit donc d’abord réécrire (2) sous la forme
(4) y10 = f1 (t, y) ; y20 = f2 (t, y) ; y1 (0) = u0 , y2 (0) = v ,
où le vecteur y a les 2 composantes y1 = u , y2 = u0 . Suivant la notation de Matlab
(voir “help ode23”), on utilise aussi t à la place de x pour la variable indépendante. En
l’occurrence, on a donc
f1 (t, y) = y2 , f2 (t, y) = c(t)y1 − f (t) .
Écrire une fonction f12(t,y) (sauvegarder dans un fichier f12.m) qui pour un scalaire
t et un vecteur y ∈ R2 calcule le vecteur des 2 composantes f1 , f2 . Pour un premier
test, on prendra l’exemple des fonctions c et f définies dans l’Exercice I.1 de la feuille
TP1. La commande
ode23(@f12,[0,1],[1,0])
affiche alors un graphe des 2 composantes de la fonction y(t) dans l’intervalle [0, 1] ,
pour les valeurs initiales y1 (0) = 1 , y2 (0) = 0 . A comparer avec la solution exacte u
donnée dans l’Exercice I.4 de la feuille TP1.
2. A partir de maintenant, on va considérer comme exemple le problème aux limites (1)
avec
c(x) = 1 ; f (x) = π 2 + 1 cos(πx) ; u0 = 1 ; u1 = 0 .
(5)
Modifier la fonction f12 pour prendre en compte ces nouvelles fonctions c et f . La
commande
[T,Y]=ode23(@f12,[0,1],[1,v]); plot(T,Y(:,1))
affiche un graphe de la solution u calculée, que l’on peut comparer avec la solution exacte
qui est de la forme
u(x) = cos(πx) + γ sinh(x)
avec une constante γ ∈ R convenable. Essayez plusieurs valeurs de v , entre autres v = 0
et v = 1 .
3. Écrire maintenant une fonction FF(v) qui utilise f12 et ode23 pour calculer une
approximation de la valeur au bord u(1) .
4. Programmer une méthode de dichotomie dicho(@F,a,b,tol) qui, à partir d’un
intervalle initial [a, b] , détermine une racine de la fonction F avec une erreur bornée par
tol. Tester avec r=dicho(@cos,0,3,1e-7), ce qui devrait donner un résultat r
satisfaisant | cos r| < 10−7 .
5. Utiliser v=dicho(@FF,0,1,1e-5) pour calculer la pente initiale v = u0 (0) cor-
respondant à la solution du problème aux limites (1) avec les données (5). Calculer cette
solution U en utilisant ode23 et dessiner le graphe de U avec la solution exacte uex .
En regardant le graphe de U − uex , noter une estimation visuelle de l’erreur maximale.
6. Adaptez la fonction errmax de la feuille TP1 à la situation actuelle et calculez
ainsi une solution du problème aux limites (1) avec les données (5) par la méthode des
différences finies qui aura la même erreur maximale que celle trouvée avec la méthode de
tir dans la question précédente. Vous avez besoin de quel nombre de nœuds N ? Compa-
rez aussi avec la méthode de Numerov (TP1/III).
7. En utilisant la fonction odeset avec les propriétés ’RelTol’ et ’AbsTol’ et
l’argument options de la fonction ode23, essayez de calculer une approximation de
la solution du problème aux limites avec une erreur inférieure à 10−9 . Comparez encore
avec les méthodes de différences finies. Connaissant le comportement asymptotique de
l’erreur des schémas de différences finies, vous pouvez extrapoler à partir d’un calcul avec
une petite valeur de N et vous n’avez ainsi pas besoin d’effectuer des calculs très gros et
longs. En comparant les temps de calcul, quelle méthode vous paraı̂t la plus efficace pour
obtenir cette précision ?