Soit \(\varepsilon=0.2\). Utiliser la méthode d’Euler explicite puis les méthodes d’Euler implicite, de Crank-Nicolson et la méthode RK4 pour calculer la solution approchée avec un pas \(h=0.1\) et comparer les deux résultats en \(t=1\). Même exercice avec \(h=0.05\) et \(h=0.01\).
Répéter avec \(\varepsilon=0.1\) puis \(\varepsilon=0\).
Commenter les résultats (hint: tracer le champ de vecteurs).
Source : (page 438 et 447) Boyce and DiPrima Elementary Differential Equations and Boundary Value Problems, seventh edition
1.1 Correction
Code
# Réinitialisation et configuration%reset -f%matplotlib inlineimport numpy as npimport matplotlib.pyplot as pltfrom scipy.optimize import fsolve# Configuration des paramètres d'affichageplt.rcParams.update({<span class="st">'font.size'</span>: <span class="dv">14</span>})# Définition de l'équation différentiellephi =lambda t,y : y**2+ t**2# Schéma d'Euler Explicite (EE)def euler_explicite(phi, tt, y0): h = tt[1] - tt[0] uu = np.zeros_like(tt) uu[0] = y0 for i inrange(len(tt) -1): uu[i +1] = uu[i] + h * phi(tt[i], uu[i]) return uu# Schéma d'Euler Implicite (EI)def euler_implicite(phi, tt, y0): h = tt[1] - tt[0] uu = np.zeros_like(tt) uu[0] = y0 for i inrange(len(tt) -1): uu[i +1] = fsolve(lambda x: -x + uu[i] + h * phi(tt[i +1], x), uu[i])[0]return uu# Schéma de Crank-Nicolson (CN)def crank_nicolson(phi, tt, y0): h = tt[1] - tt[0] uu = np.zeros_like(tt) uu[0] = y0 for i inrange(len(tt) -1): f =lambda x: -x + uu[i] + h * (phi(tt[i], uu[i]) + phi(tt[i +1], x)) /2 uu[i +1] = fsolve(f, uu[i])[0]return uu# Schéma de Runge-Kutta d'ordre 4 (RK4)def runge_kutta_4(phi, tt, y0): h = tt[1] - tt[0] uu = np.zeros_like(tt) uu[0] = y0 for i inrange(len(tt) -1): k1 = phi(tt[i], uu[i]) k2 = phi(tt[i] + h /2, uu[i] + h * k1 /2) k3 = phi(tt[i] + h /2, uu[i] + h * k2 /2) k4 = phi(tt[i +1], uu[i] + h * k3) uu[i +1] = uu[i] + (k1 +2* k2 +2* k3 + k4) * h /6return uu# Initialisation des paramètrest0, tfinal, y0 =0, 1, 1# Création d'un tableau pour stocker tous les résultatsresults = []# Création de la figure avec 3 lignes et 3 colonnesfig, axes = plt.subplots(3, 3, figsize=(18, 18)) for row, epsilon inenumerate([0.2, 0.1, 0]):for col, h inenumerate([0.1, 0.05, 0.01]): N =int((1- epsilon) / h) tt = np.linspace(0, 1- epsilon, N +1) # Utilisation de linspace pour éviter les erreurs d'arrondi uu_EE = euler_explicite(phi, tt, y0) uu_EI = euler_implicite(phi, tt, y0) uu_CN = crank_nicolson(phi, tt, y0) uu_RK4 = runge_kutta_4(phi, tt, y0)# Enregistrer les résultats dans le tableau results.append([epsilon, h, uu_EE[-1], uu_EI[-1], uu_CN[-1], uu_RK4[-1]])# Tracer les solutions dans le sous-graphe correspondantif row ==0: axes[row, col].set_title(f"$h={</span>h<span class="sc">}$", fontsize=14) axes[row, col].plot(tt, uu_EE, '*-', label=f'EE', color='r') axes[row, col].plot(tt, uu_EI, 'o-', label=f'EI', color='g') axes[row, col].plot(tt, uu_CN, 's-', label=f'CN', color='m') axes[row, col].plot(tt, uu_RK4, 'D-', label=f'RK4', color='b') axes[row, col].set_xlabel('$t$', fontsize=12) axes[row, col].set_ylabel('$y(t)$', fontsize=12) axes[row, col].set_xlim([0, 1]) axes[row, col].grid(True) axes[row, col].legend(loc='best') axes[row, 0].set_ylabel(f'tfinal=$={</span><span class="dv">1</span> <span class="op">-</span> epsilon<span class="sc">}$', fontsize=16) # Ajouter l'étiquette de epsilonplt.tight_layout() # Pour un espacement optimal entre les sous-figuresplt.show()# Afficher tous les résultats dans un tableaufrom tabulate import tabulateheaders = ["tfinal", "h", "EE(tfinal)", "EI(tfinal)", "CN(tfinal)", "RK4(tfinal)"]table = []# Remplir le tableau avec les résultats, en ajoutant un séparateur après chaque changement de tfinalfor result in results: table.append([1- result[0], result[1], result[2], result[3], result[4], result[5]])table_str = tabulate(table, headers=headers, tablefmt="html", floatfmt=".2f")display(table_str)
/tmp/ipykernel_84147/1899722127.py:13: RuntimeWarning: overflow encountered in scalar power
phi = lambda t,y : y**2 + t**2
tfinal
h
EE(tfinal)
EI(tfinal)
CN(tfinal)
RK4(tfinal)
0.80
0.10
3.51
5.00
7.45
5.84
0.80
0.05
4.20
10.00
6.08
5.85
0.80
0.01
5.34
6.57
5.86
5.85
0.90
0.10
4.80
5.00
10.00
14.02
0.90
0.05
6.46
10.00
20.01
14.27
0.90
0.01
10.80
28.12
14.44
14.30
1.00
0.10
7.19
5.00
10.00
735.10
1.00
0.05
12.32
10.00
20.00
175862.67
1.00
0.01
90.76
50.00
100.00
inf
Le problème est que nous avons présupposé que le problème est mathématiquement bien posé, c’est-à-dire qu’une unique solution existe sur l’intervalle d’étude \([0;1]\). Il se trouve qu’aucune solution n’existe sur l’interval \([0;1]\).
La théorie générale permet de montrer qu’une solution unique existe pour un certain intervalle contenant \(t=0\), mais elle ne nous dit pas jusqu’où s’étend cet intervalle. Nous pouvons montrer qu’une solution existe pour au moins \(t>\pi/4\). Cependant, la solution devient non bornée quelque part entre \(\pi/4\) et \(1\). Nous pouvons le prouver en observant que la solution cherchée \(y\) est bornée par les deux fonctions \(z\) et \(w\) suivantes
\(z(t)=\frac{1}{1-t}\) qui est définie sur \([0;1[\),
\(w(t)=\tan\left(t+\frac{\pi}{4}\right)\) qui est définie sur \([0;\pi/4[\).
Puisque \(z(t)\to+\infty\) quand \(t\to 1^-\), la solution cherchée \(y\) devient non bornée quelque part entre \(\pi/4\) et \(1\).
Même sans connaître la solution exacte, nous pouvons étudier le comportement sur l’intervalle \([0;1]\) en affichant le champ de vecteurs :
Code
figure(figsize=(10,7))# quiverplot# define a grid and compute direction at each pointg1 = linspace(0,1,21)g2 = linspace(0,3,21)T,Y = meshgrid(g1,g2) # create a gridDT, DY =1, phi(T,Y) # compute growth rate on the gridM = sqrt(DT**2+DY**2) # norm growth rate M[ M==0 ] =1# avoid zero division errors quiver(T,Y, DT/M, DY/M, M, pivot='mid')# streamplot(T,Y, DT/M, DY/M, color=M)grid()xlabel('t')ylabel('y')title('Champ des pentes');
Cet exemple montre l’importance de l’interaction entre la théorie et le calcul numérique.
Si vous vous contentez de calculer numériquement une solution approchée avec un schéma et un seul pas \(h\) sans savoir qu’une solution unique existe, vous pouvez penser que vous avez réussi alors que vous êtes en train de calculer quelque chose qui n’existe pas.
Lorsque cela est possible, il faut essayer d’abord de prouver que le problème est mathématiquement bien posé avant d’en calculer une solution approchée. Par exemple, ce serait bien de connaître la taille de l’intervalle d’existence de la solution avant de calculer quoi que ce soit. Bien-sûr, ce n’est pas toujours possible. Dans notre exemple il n’est pas évident, en regardant l’équation, de voir qu’il y a un problème pour \(t \to 1\). L’affichage du champ de vecteurs peut aider à voir ce genre de problème.
Une bonne pratique numérique, c’est-à-dire l’essai de plus d’une taille de pas et/ou de plus d’une méthode numérique, est également précieuse. Les difficultés que nous avons rencontrées avec les comparaisons suggèrent qu’il pourrait y avoir un problème théorique sous-jacent qu’il faut clarifier.
On ne peut pas dire que la théorie ou le calcul viennent nécessairement en premier, il faut itérer entre les deux, en commençant par l’approche la plus facile à mettre en œuvre. Les résultats théoriques sont plus satisfaisants lorsqu’ils sont disponibles, mais la théorie ne nous en dit souvent pas autant que ce que nous aimerions savoir. En outre, les gens font des erreurs dans les calculs théoriques, tout comme dans les calculs numériques. Il est préférable que les travaux théoriques et numériques se valident mutuellement.
Ce problème montre l’importance de se préoccuper de l’existence et de l’unicité, mais les méthodes théoriques ne sont pas les seules méthodes pour explorer l’existence. Quoi qu’il en soit, le problème montre que sans une certaine diligence - théorique ou numérique - vous pourriez naïvement calculer une “solution” approximative là où il n’existe aucune solution.
2 Exercice : Étude d’un problème numériquement mal posé
Utiliser la méthode RK4 pour calculer la solution approchée avec un pas \(h=\frac{1}{10}\) et comparer à la solution exacte pour cinq valeurs des \(\varepsilon\), à savoir \(\varepsilon=\pm 1\), \(\varepsilon=\pm \frac{1}{10}\) et \(\varepsilon=0\).
Même exercice avec la méthode d’Euler explicite.
Commenter les résultats.
2.1 Correction : solution exacte
La solution exacte est \(y(t)=\varepsilon e^{3t}+e^{-t}\) (ci-dessous, le calcul de la solution exacte avec sympy).
Si \(\varepsilon=0\), la solution exacte devient donc \(y(t)=e^{-t}\) mais le problème est mal conditionné : toute (petite) erreur de calcul a le même effet qu’une perturbation de la condition initiale. En d’autres termes, on “réveil” le terme dormant \(e^{3t}\).
Affichons l’allure de la solution exacte dans les cinq cas \(\varepsilon=\pm1\), \(\varepsilon=\pm\frac{1}{10}\) et \(\varepsilon=0\). On remarque que la solution exacte est très sensible à la condition initiale.
3 Exercice : A-stabilité du schéma d’Euler expicite dans le cas d’un système
Soit le problème de Cauchy:
\(\begin{cases}
y''(t)+1001y'(t)+1000y(t)=0, & \forall t \in \mathbb{R},\\
y(0)=1,\\
y'(0)=0.
\end{cases}\)
Écrire le problème comme un système de 2 EDO d’ordre 1 et en calculer la solution.
Écrire le schéma le schéma d’Euler explicite pour ce problème.
Pour quelles valeurs de \(h\) le schéma est A-stable?
3.1 Correction
Soit \(\mathbf{z}(t)=(y(t),y'(t))^T\) alors \(y''(t)+1001y'(t)+1000y(t)=0\) se réécrit \(
\mathbf{z}'(t)=
-\begin{pmatrix}
0&-1\\
1000&1001
\end{pmatrix}
\mathbf{z}(t)
\) Cette matrice a pour valeurs propres \(\lambda_1=1\) et \(\lambda_2=1000\) et pour vecteurs propres associés \(\mathbf{v}_1=(-1,1)^T\) et \(\mathbf{v}_2=(0,1000)^T\) (ci dessous le calcul formel des valeurs et vecteurs propes). Ainsi \(
\mathbf{z}(t)=c_1\mathbf{v}_1e^{-\lambda_1t}+c_2\mathbf{v}_2e^{-\lambda_2t}
\) et \(
y(t)=-c_1e^{-t}-c_2e^{-1001t}.
\) Il s’agit d’un problème stiff!
Code
%reset -fimport sympy as symbsymb.init_printing()symb.var('t')M = symb.Matrix(((0,-1), (1000,1001)))display(M)p=M.charpoly(t)display(symb.factor(p))P, D = M.diagonalize()display(P,D)display(symb.simplify(P*D*P**-1))
Le schéma est A-stable si \(\lim_{n\to+\infty}\mathbf{u}_{n+1}=\mathbf{0}\). Il faut donc que \(\begin{cases}
|1-1000h|<1\\
|1-h|<1
\end{cases}\) d’où la condition \(h<\frac{2}{1000}\).
5 Exercice : Étude théorique et numérique du mouvement d’un pendule
Considérons le problème du mouvement d’un pendule de masse \(m\) suspendu au point \(O\) par un fil non pesant de longueur \(\ell\).
On se propose d’étudier l’évolution au cours du temps \(t\ge0\) de l’angle \(\vartheta(t)\) autour du point \(0\).
L’angle \(\vartheta(t)\) est mesuré par rapport à une verticale passante par \(0\).
Considérons pour simplifier seulement la force de gravité \(g\).
La fonction \(\vartheta\) est alors solution de l’EDO d’ordre 2
Montre analytiquement que l’énergie totale reste constante au cours du temps, i.e. que \(\frac{\mathrm{d}}{\mathrm{d}t}E(q(t),p(t))=0\).
Tracer avec Python les courbes de niveau de la fonction \((q,p)\mapsto E(q,p)\), i.e. les ensembles \(\mathcal{N}_\kappa = \{(q,p)\in[-\pi;\pi]\times\mathbb{R}\ |\ E(q,p)=\kappa\} =\{(q,p)\in[-\pi;\pi]\times\mathbb{R}\ |\ p^2=2(\kappa+\cos(q))\}. \)
Lorsque \(\vartheta\) est petit on peut considérer l’approximation \(\sin(\vartheta)\simeq\vartheta\). Calculer la solution exacte de cette équation approchée.
Étude numérique
Dans une simulation numérique on aimerait que la propriété de conservation de l’énergie soit préservée aussi bien que possible. Nous allons regarder ce qui se passe avec certaines méthodes vues en cours.
On notera \(x_n\approx q(t_n)=\vartheta(t_n)\) et \(y_n\approx p(t_n)=\vartheta'(t_n)\).
À chaque instant \(t_n\), on calculera les valeurs approchées de l’énergie cinétique \(E_c(t_n)=\frac{1}{2}y_n^2\), de l’énergie potentielle \(E_p(t_n)=-\cos(x_n)\) et de l’énergie totale \(E(t_n)=E_c(t_n)+E_p(t_n)\).
Choisissons \(h=t_{n+1}-t_n=0.1\) et \((x_0,y_0)=(\pi/4,0)\).
Dans un premier temps on se propose d’appliquer la méthode d’Euler explicite à la résolution du système.
Tracer \((t_n,x_n)\) et \((t_n,y_n)\) sur un même graphique;
sur un autre graphique représenter les trois énergies simultanément en fonction du temps;
tracer enfin \((x_n,y_n)\) sur un autre graphique et vérifier que la solution numérique tourne vers l’extérieur (i.e. l’énergie augmente au cours du temps).
Considérons le schéma obtenu en écrivant le schéma d’Euler explicite pour la première équation et le schéma d’Euler implicite pour la seconde.
Montrer que ce schéma peut s’écrire sous forme parfaitement explicite;
tracer \((t_n,x_n)\) et \((t_n,y_n)\) sur un même graphique;
sur un autre graphique représenter les trois énergies simultanément en fonction du temps;
tracer enfin \((x_n,y_n)\) sur un autre graphique. Que peut-on constater? Commenter les résultats.
On se propose d’appliquer les méthodes d’Euler modifiée, de Heun, AB2, AB3 et RK4 à la résolution du système. Pour chaque méthode:
Tracer \((t_n,x_n)\) et \((t_n,y_n)\) sur un même graphique;
sur un autre graphique représenter les trois énergies simultanément en fonction du temps;
tracer enfin \((x_n,y_n)\) sur un autre graphique. Que peut-on constater? Commenter les résultats.
5.1 Correction étude théorique
Il s’agit d’un système de deux EDO scalaires d’ordre 1 du type
L’énergie totale reste constante au cours du temps: \(\begin{aligned}
\frac{\mathrm{d}}{\mathrm{d}t}E(q(t),p(t))
&= \sin(q(t))q'(t)+p(t)p'(t)
{}= \sin(q(t))p(t)+p(t)p'(t) \\
&= p(t)\left( \sin(q(t))+p'(t) \right)
{}= p(t)\left( \sin(q(t))-\sin(q(t)) \right)=0.
\end{aligned}\)
Courbes de niveau de la fonction \((q,p)\mapsto E(q,p)\):
\(\vartheta''(t)+\omega^2\vartheta(t)=0\) est une EDO d’ordre 2 homogène à coefficients constants dont la solution est \(\vartheta(t)=\vartheta(0)\cos(\omega t)+\frac{\vartheta'(0)}{\omega}\sin(\omega t)\).
def euler_progressif(phi1,phi2,tt,y0): h = tt[1]-tt[0] q = [y0[0]] p = [y0[1]]for i inrange(len(tt)-1): q.append(q[i]+h*phi1(tt[i],q[i],p[i])) p.append(p[i]+h*phi2(tt[i],q[i],p[i]))return [q,p][q_ep, p_ep] = euler_progressif(phi1,phi2,tt,y0)Ec_ep = [p**2/2for p in p_ep]Ep_ep = [-math.cos(q) for q in q_ep]Et_ep = [Ec_ep[i]+Ep_ep[i] for i inrange(len(Ec_ep))]print ('Euler explicite : |energie totale (t=Tfinal) - Energie totale (t=0)|=', abs(Et_ep[-1]-Et_ep[0]))figure(figsize=(18,5))subplot(131)plot(tt,q_ep,tt,p_ep)xlabel('t')legend(['x(t)','y(t)'])title('Euler progressif - x(t) et y(t)')subplot(132)plot(q_ep,p_ep)xlabel('x(t)')ylabel('y(t)')title('Euler progressif')subplot(133)plot(tt,Ec_ep,tt,Ep_ep,tt,Et_ep)xlabel('t')legend(['Cinetique','Potentielle','Totale'])title('Euler progressif - Energies');
Euler explicite : |energie totale (t=Tfinal) - Energie totale (t=0)|= 3.1764943624320936
5.2.2 Méthode d’Euler explicite/implicite
Schéma obtenu en écrivant le schéma d’Euler explicite pour la première équation et le schéma d’Euler implicite pour la seconde:
Dans notre cas \(\varphi_2(t_{n+1},x_{n+1},y_{n+1})=-\sin(x_{n+1})\) donc le schéma est parfaitement explicite:
Code
def euler_symplectique(phi1,phi2,tt,y0): h = tt[1]-tt[0] q = [y0[0]] p = [y0[1]]for i inrange(len(tt)-1): q.append(q[i]+h*phi1(tt[i],q[i],p[i])) p.append(p[i]+h*phi2(tt[i+1],q[i+1],p[i])) # il faudrait p[i+1] (schema implicite) mais comme dans notre cas phi2 ne depend pas de "p", le schema est explicitereturn [q,p][q_es, p_es] = euler_symplectique(phi1,phi2,tt,y0)Ec_es = [p**2/2for p in p_es]Ep_es = [-math.cos(q) for q in q_es]Et_es = [Ec_es[i]+Ep_es[i] for i inrange(len(Ec_es))]print ('Euler symplectique : |energie totale (t=Tfinal)-Energie totale (t=0)|=', abs(Et_es[-1]-Et_es[0]))figure(figsize=(18,5))subplot(131)plot(tt,q_es,tt,p_es)xlabel('t')legend(['x(t)','y(t)'])title('Euler symplectique - x(t) et y(t)')subplot(132)plot(q_es,p_es)xlabel('x(t)')ylabel('y(t)')title('Euler symplectique - y(x)')subplot(133)plot(tt,Ec_es,tt,Ep_es,tt,Et_es)xlabel('t')legend(['Cinetique','Potentielle','Totale'])title('Euler symplectique - Energies');
Euler symplectique : |energie totale (t=Tfinal)-Energie totale (t=0)|= 0.030048843024607974
