TD1 - Exercices sur les schémas à un pas

Exercice : Modèle de Malthus

On modélise la croissance d’une population d’effectif \(y(t)\ge0\) par l’équation différentielle

\(\begin{equation} y'(t)=a y(t), \qquad a\in\mathbb{R}. \end{equation}\)

Dans ce modèle la variation de population est proportionnelle à la population elle-même : le coefficient \(a\) est appelé paramètre de Malthus ou potentiel biologique et représente la différence entre le nombre de nouveaux nés par individu en une unité de temps et la fraction d’individus qui meurent en une unité de temps. On s’intéresse à l’unique solution de cette équation vérifiant \(y(t_0)=y_0\) où \(y_0\) et \(t_0\) sont des réels donnés.

1. Donner l’unique solution théorique à ce problème de Cauchy. Calculer \(\lim_{t\to+\infty}y(t)\) en fonction du paramètre \(a\).
2. Calculer analytiquement la solution obtenue avec la méthode d’Euler explicite sur l’intervalle \([0;2]\) avec \(N_h=5\), ensuite la tracer avec Python et la comparer à la solution exacte.
3. Soit \(u_{N_h}(2)\) l’approximation de \(y(2)\) obtenue avec la méthode d’Euler explicite et \(N_h\) nœuds. Calculer \(u_{N_h}(2)\) avec \(N_h=2^i\) et montrer analytiquement que \(\lim_{N_h\to+\infty}u_{N_h}(2)=y(2)\).
4. Tracer le graphe de \(\max_{0\le n\le N_h}|y(t_n)-u_n|\) en fonction de \(N_h\) ainsi que le graphe des courbes \(1/N_h\), \(1/N_h^2\) … avec une progression logarithmique en abscisse et en ordonnées. On prendra par exemple \(N_h=200\), puis \(N_h=300\), \(400\), … \(2000\).

Correction

L’unique solution du problème de Cauchy \(\begin{cases}y'(t)=a y(t)\\y(t_0)=y_0>0\end{cases}\) pour \(t\ge t_0\) est la fonction \(y(t)=y_0e^{at}\). On a

\( \lim_{t\to+\infty}y(t)= \begin{cases} +\infty &\text{si } a>0,\\ y_0 &\text{si } a=0,\\ 0^+ &\text{si } a<0, \end{cases} \)

Pour \(a=1\), \(t_0=0\) et \(y_0=1\), la solution exacte est \(y(t)=e^t\).

On subdivise l’intervalle \([t_0;T]\) en \(N\) intervalles \([t_{n};t_{n+1}]\) de largeur \(h=(T-t_0)/N\) avec \(t_n=t_0+nh\) pour \(n=0,\dots,N\).
Pour chaque nœud \(t_n\), on note \(y_n=y(t_n)\); l’ensemble des valeurs \(\{y_0, y_1,\dots , y_{N} \}\) représente la solution exacte discrète.
Pour chaque nœud \(t_n\), on cherche la valeur inconnue \(u_n\) qui approche la valeur exacte \(y_n\).
L’ensemble des valeurs \(\{u_0 = y_0, u_1,\dots , u_{N} \}\) représente la solution numérique.

Dans la méthode d’Euler explicite, cette solution approchée est obtenue en construisant une suite récurrente comme suit

\(\begin{cases} u_0=y_0,\\ u_{n+1}=u_{n}+h\varphi(t_{n},u_{n}). \end{cases}\)

En utilisant ce schéma pour notre EDO on obtient

\(\begin{cases} u_0=y_0,\\ u_{n+1}=(1+ah)u_n \end{cases}\)

soit encore \(u_n=(1+ah)^{n}u_0=\left(1+\frac{2}{5}\right)^n\).

import matplotlib.pyplot as plt
import numpy as np
plt.rcParams.update({'font.size': 16})


def EE(phi, tt, y0):
    uu = np.empty_like(tt)
    uu[0] = y0
    h = tt[1] - tt[0]
    for i in range(len(tt) - 1):
        uu[i + 1] = uu[i] + h * phi(tt[i], uu[i])
    return uu

phi = lambda t, y: a * y  
sol_exacte = lambda t: y0 * np.exp(a * t)

# INITIALISATION
a = 1
N = 5
t0, y0 = 0, 1
tfinal = 2

# CALCUL SOLUTION APPROCHEE
tt = np.linspace(t0, tfinal, N + 1)
uu = EE(phi, tt, y0)

# CALCUL SOLUTION EXACTE POUR PLOT
yy = sol_exacte(tt)   # [sol_exacte(t) for t in tt]

# PLOT
plt.plot(tt, yy, 'm-', label='Exacte')
plt.plot(tt, uu, 'go-', label='Euler explicite')
plt.legend()
plt.show()

# TABLEAU
from tabulate import tabulate
T = []
T.append(["n", "y_n", "u_n", "(1+h)^n"])
h1 = tt[1] - tt[0]
for n in range(N):
    T.append([n, yy[n], uu[n], (1 + h1) ** n])
display(tabulate(T, headers="firstrow", tablefmt="html"))

n	y_n	u_n	(1+h)^n
0	1	1	1
1	1.49182	1.4	1.4
2	2.22554	1.96	1.96
3	3.32012	2.744	2.744
4	4.95303	3.8416	3.8416

On compare les approximations de \(y(2)\) avec plusieurs valeurs de \(N_h\) :

\( u_{N_h}=(1+ah)^{N_h}u_0=\left(1+\frac{2}{N_h}\right)^{N_h}\xrightarrow[N_h\to+\infty]{}e^2=y(2). \)

from IPython.display import Markdown

y_2 = sol_exacte(tfinal)
display(Markdown(f'Solution exacte en $t=2$ : $y(2) = {y_2}$'))
display(Markdown('Solution approchee en $t=2$ : $u(2)$ avec $N$ points:'))
from tabulate import tabulate
table_data = []
table_data.append(["N", "u(2)", "|y(2)-u(2)|"])
for i in range(3, 15):
    N = 2 ** i
    tt = np.linspace(t0, tfinal, N + 1)
    uu = EE(phi, tt, y0)
    error = abs(uu[-1] - y_2)
    table_data.append([N, f'{uu[-1]:1.8f}', f'{error:1.8f}'])
display(tabulate(table_data, headers="firstrow", tablefmt="html"))

Solution exacte en \(t=2\) : \(y(2) = 7.38905609893065\)

Solution approchee en \(t=2\) : \(u(2)\) avec \(N\) points:

N	u(2)	|y(2)-u(2)|
8	5.96046	1.42859
16	6.58325	0.805806
32	6.95867	0.430389
64	7.16628	0.22278
128	7.27567	0.113386
256	7.33185	0.0572055
512	7.36032	0.0287325
1024	7.37466	0.0143989
2048	7.38185	0.00720766
4096	7.38545	0.00360588
8192	7.38725	0.00180346
16384	7.38815	0.00090186

Vérifions la convergence linéaire:

error = []
H = []
N = np.arange(200, 1000, 100)
for n in N:
    tt = np.linspace(t0, tfinal, n + 1)
    H.append(tt[1] - tt[0])
    yy = sol_exacte(tt)
    uu = EE(phi, tt, y0)
    error.append(np.linalg.norm(uu-yy, np.inf))

# ESTIMATION ORDRES DE CONVERGENCE par regression
pente = np.polyfit(np.log(H), np.log(error), 1)[0]

# PLOT ERREUR  
plt.figure(figsize=(10, 6))
plt.title(f'Ordre = pente obtenue par régression linéaire = {pente}')
plt.loglog(H, error, '*-');

Exercice : Modèle de Verhulst (croissance logistique).

Dans le modèle logistique, on suppose que la croissance de la population d’effectif \(y(t)\ge0\) est modélisée par l’équation différentielle :

\(\begin{equation} y'(t)=(a-by(t))y(t), \qquad a,b>0. \end{equation}\)

Ce modèle prend en compte une éventuelle surpopulation. En effet, le taux de croissance de la population est désormais de la forme \(a-by(t)\): c’est la différence entre un taux de fertilité \(a\) et un taux de mortalité de la forme \(by(t)\) qui augmente avec l’effectif de la population.

1. Calculer les solutions exactes de l’équation différentielle.

2. Écrire la relation entre \(u_{n+1}\) et \(u_n\) obtenue en appliquant la méthode d’Euler explicite.

3. Pour \(a=b=0.5\) et \(t_0=0\), tracer avec python les solutions exactes sur l’intervalle de temps \([0;20]\). Sur un autre graphique à côté, tracer les solutions correspondantes obtenues par la méthode d’Euler explicite avec \(N_h=50\). Sur un autre graphique à côté, tracer les solutions correspondantes obtenues par la méthode d’Euler explicite avec \(N_h=100\). Considérer les valeurs de \(y_0\) suivantes :

   \(\begin{aligned} y_0&=0.01 & y_0&=0.1 & y_0&=0.4 & y_0&=0.8 & y_0&=1 & y_0&=1.5 & y_0=2 \end{aligned} \)

Correction

On cherche l’unique solution du problème de Cauchy \(\begin{cases}y'(t)=(a-by(t))y(t)\\y(t_0)=y_0\ge0\end{cases}\) pour \(t\ge t_0\).

On écrit l’EDO sous la forme \(y'(t) = -b y(t) \big( y(t)-\frac{a}{b} \big).\)

Si \(y(t)=A\) pour tout \(t>0\) est solution de l’EDO, alors on doit avoir \(A=(a-bA)A\). Donc, soit \(A=0\) soit \(A=\frac{a}{b}\).
Ainsi la solution constante \(y(t)=\frac{a}{b}>0\) représente l’effectif d’équilibre de la population, directement liée aux contraintes environnementales (densité, ressources, etc.).

De plus, les solutions sont des fonctions strictement croissantes si \(y(t)\in]0,a/b[\) et décroissantes si \(y(t)>a/b\).

Calculons la solution exacte en remarquant qu’il s’agit d’une EDO à variables séparables. Pour \(y(t)\neq0\) et \(y(t)\neq\frac{a}{b}\), on remarque que \( -b=\frac{y'(t)}{ \left( y(t)-\frac{a}{b} \right) y(t) } =\frac{\frac{b}{a}y'(t)}{y(t)-\frac{a}{b}y'(t)}-\frac{\frac{b}{a}}{y(t)}. \) On doit alors calculer \( \int\frac{1}{y-\frac{a}{b}}\mathrm{d}y-\int\frac{1}{y}\mathrm{d}y=-a\int 1 \mathrm{d}t. \) On obtient \( \ln\left( 1-\frac{b/a}{y} \right)=-at+C \) et on en déduit \( y(t)=\frac{a/b}{1-De^{-at}}. \) En imposant la condition initiale \(y(0)=y_0\) on trouve la constante d’intégration \(D\) et on conclut que toutes les solutions du problème de Cauchy pour \(t\ge0\) sont

\( y(t) = \begin{cases} 0 & \text{si } y_0=0,\\(0.5em] \dfrac{a}{b} & \text{si } y_0=\dfrac{a}{b},\\(0.5em] \dfrac{1}{\left(\frac{1}{y_0}-\frac{b}{a}\right) e^{-at}+\frac{b}{a}} &\text{sinon.} \end{cases} \)

Remarquons que \(\lim_{t\to+\infty}y(t)=\frac{a}{b}\): une population qui évolue à partir de \(y_0\) individus à l’instant initiale tend à se stabiliser vers un nombre d’individus d’environ \(a/b\), ce qui représente la capacité de l’environnement.

La méthode d’Euler explicite (ou progressive) est une méthode d’intégration numérique d’EDO du premier ordre de la forme \(y'(t)=\varphi(t,y(t))\).
On subdivise l’intervalle \([t_0;T]\) en \(N\) intervalles \([t_{n};t_{n+1}]\) de largeur \(h=\dfrac{T-t_0}{N}\) avec \(t_n=t_0+nh\) pour \(n=0,\dots,N\). La longueur \(h\) est appelé le pas de discrétisation.
Pour chaque nœud \(t_n\), on note \(y_n=y(t_n)\); l’ensemble des valeurs \(\{y_0, y_1,\dots , y_{N} \}\) représente la solution exacte discrète.
Pour chaque nœud \(t_n\), on cherche la valeur inconnue \(u_n\) qui approche la valeur exacte \(y_n\).
L’ensemble des valeurs \(\{u_0 = y_0, u_1,\dots , u_{N_h} \}\) représente la solution numérique.

Dans la méthode d’Euler explicite, cette solution approchée est obtenue en construisant une suite récurrente comme suit \(\begin{cases} u_0=y_0,\\ u_{n+1}=u_{n}+h\varphi(t_{n},u_{n}). \end{cases}\) En utilisant cette construction pour notre EDO on obtient \(\begin{cases} u_0=y_0,\\ u_{n+1}=u_{n}+h(a-bu_n)u_n \end{cases}\)

import numpy as np
import matplotlib.pyplot as plt

# Réinitialisation des paramètres par défaut de Matplotlib
plt.rcdefaults()

# Mise à jour des paramètres de police
plt.rcParams.update({'font.size': 12})
plt.rcParams.update({'figure.max_open_warning': 0})


def EE(phi, tt, y0):
    h = tt[1] - tt[0]
    uu = np.empty_like(tt)
    uu[0] = y0
    for i in range(len(tt) - 1):
        uu[i + 1] = uu[i] + h * phi(tt[i], uu[i])
    return uu

phi = lambda t, y: (a - b * y) * y


def sol_exacte(t, y0):
    if abs(y0 - 0) <= 1.e-10:  # avec des float on ne peut pas utiliser ==0
        return 0
    elif abs(y0 - a / b) <= 1.e-10:
        return a / b
    else:
        return 1 / ((1 / y0 - b / a) * np.exp(-a * t) + b / a)

# Vectoriser la fonction sol_exacte
sol_exacte_vectorized = np.vectorize(sol_exacte, excluded=['y0'])


# INITIALISATION
a = 0.5
b = 0.5
t0 = 0
tfinal = 20
yy_0 = [0.001, 0.1, 0.5, 1, 1.5, 2]
NN = [25, 50, 100]


# MAIN
f, axs = plt.subplots(1, 3, figsize=(15, 5))
axs = np.ravel(axs)
colors = plt.cm.Set1(np.linspace(0, 1, len(yy_0)))
for ax,N in zip(axs,NN):
    tt = np.linspace(t0, tfinal, N)
    for i, y0 in enumerate(yy_0):
        color = colors[i]
        uu = EE(phi, tt, y0)
        ax.plot(tt, uu, 'o', label=f'Approchée avec y_0={y0}',color=color)
        yy = sol_exacte_vectorized(tt, y0)
        ax.plot(tt, yy, label=f'Exacte avec y_0={y0}',color=color)
    ax.set_title(f'{N} noeuds')
    ax.set_xlim(t0, tfinal)
    ax.set_ylim(0, 2.)
    ax.grid()
    # ax.legend()

plt.suptitle('En continue les solutions exactes et en ronds les solutions approchées')
plt.tight_layout()
plt.show()

Exercice : concentration

L’évolution de la concentration de certaines réactions chimiques au cours du temps peut être décrite par l’équation différentielle

\( y'(t)=-\frac{1}{1+t^2}y(t). \)

Sachant qu’à l’instant \(t=0\) la concentration est \(y(0)=5\), déterminer la concentration à \(t=2\) à l’aide de la méthode d’Euler implicite avec un pas \(h=0.5\).

Correction

Dans ce cas, le schéma d’Euler implicite s’écrit

\( u_{n+1}=u_n-\frac{h}{1+t_{n+1}^2}u_{n+1}. \)

Elle peut être résolue algébriquement et cela donne la suite

\( u_{n+1}=\frac{u_n}{1+\frac{h}{1+t_{n+1}^2}}. \)

Si à l’instant \(t=0\) la concentration est \(y(0)=5\), et si \(h=1/2\), alors \(t_n=n/2\) et

\( u_{n+1}=\frac{4+(n+1)^2}{6+(n+1)^2}u_n. \)

On obtient donc

\( \begin{array}{ccc} n & t_n & u_n \\ \hline 0 & 0 & 5 \\ 1 & 0.5 & \frac{4+1^2}{6+1^2}5 =\frac{5}{7}5 =\frac{25}{7}\approx3.57\\ 2 & 1.0 & \frac{4+2^2}{6+2^2}\frac{25}{7} =\frac{8}{10}\frac{25}{7}=\frac{20}{7}\approx2.86\\ 3 & 1.5 & \frac{4+3^2}{6+3^2}\frac{20}{7} =\frac{13}{15}\frac{20}{7}=\frac{52}{21}\approx2.48\\ 4 & 2.0 & \frac{4+4^2}{6+4^2}\frac{52}{21}=\frac{20}{22}\frac{52}{21}=\frac{520}{231}\approx2.25\\ % 5 & 2.5 & \sfrac{15080}{7161}\approx2.11\\ % 6 & 3.0 & \sfrac{301600}{150381}\approx2.01 \end{array} \)

La concentration à \(t=2\) est d’environ \(2.25\) qu’on peut comparer avec le calcul exact \(y(2)=5e^{-\arctan(2)}\approx1.652499838\).

Exercice : A-stabilité du schéma de Crank-Nicolson

Soit le problème de Cauchy:

\(\begin{cases} y'(t)+10y(t)=0, & \forall t \in \mathbb{R},\\ y(0)=y_0>0. \end{cases}\)

1. Montrer qu’il existe une unique solution globale \(y\) \(\in\) \(\mathscr{C}^\infty(\mathbb{R},\mathbb{R})\) que vous préciserez explicitement.
2. Soit \((u_n)_{n\in \mathbb{N}}\) la suite obtenue avec le schéma numérique de Crank-Nicolson. Montrer que \((u_n)_{n\in \mathbb{N}}\) est une suite géométrique dont on précisera la raison. Donner l’expression de \(u_n\) en fonction de \(n\).
3. Montrer que la raison \(r\) de la suite vérifie \(|r|<1\) pour tout \(h>0\). Ce schéma est-il inconditionnellement A-stable?
4. Sous quelle condition sur \(h>0\) le schéma génère-t-il une suite positive ?

Correction

C’est un problème deCauchy du type \(\begin{cases} y'(t)=\varphi(t,y(t)), & \forall t \in \mathbb{R},\\ y(0)=y_0>0, \end{cases}\) avec \(\varphi(t,y)=g(y)=-10y\).

On montre d’abord qu’il existe une et une seule solution locale (ie sur \([-T;T]\)) de classe \(\mathscr{C}^1([-T,T],\mathbb{R})\).
On montre ensuite que cette solution est de classe \(\mathscr{C}^\infty([-T,T],\mathbb{R})\).
On montre enfin que la solution admet un prolongement sur \(\mathbb{R}\). + Comme \(g\) \(\in\) \(\mathscr{C}^1(\mathbb{R},\mathbb{R})\), d’après le théorème de Cauchy-Lipschitz il existe \(T>0\) et une unique solution \(y\) \(\in\) \(\mathscr{C}^1([-T,T],\mathbb{R})\). + Par récurrence, en exploitant l’EDO et la régularité de \(g\), on grimpe en régularité sur \(y\) ainsi \(y\) \(\in\) \(\mathscr{C}^\infty([-T,T],\mathbb{R})\). + La fonctionne nulle est solution de l’EDO (mais non du problème de Cauchy donné). Par l’unicité de la solution du problème de Cauchy on en déduit que soit \(y(t)>0\) pour tout \(t\) \(\in\) \([-T,T]\) (ie lorsque \(y_0>0\)) soit \(y(t)<0\) pour tout \(t\) \(\in\) \([-T,T]\) (ie lorsque \(y_0>0\)). De plus, \(y\) est décroissante si \(y_0>0\) et croissante si \(y_0<0\). On en déduit par le théorème des extrémités que la solution \(u\) admet un prolongement sur \(\mathbb{R}\) solution de l’EDO.

Pour en calculer la solution, on remarque qu’il s’agit d’une EDO à variables séparables. L’unique solution constante est \(y(t)\equiv0\), toutes les autres solutions sont du type \(y(t)=C e^{-10t}\). En prenant en compte la condition initiale on conclut que l’unique solution du problème de Cauchy est \( y(t)=y_0e^{-10t} \) définie pour tout \(t\) \(\in\) \(\mathbb{R}\).

import matplotlib.pyplot as plt
import numpy as np

# Réglage de la taille de la police
plt.rcParams.update({'font.size': 12})

# Définition de la fonction exacte
y0 = 1
exacte = lambda t: y0 * np.exp(-10*t)

# Génération des données pour le tracé
tt = np.linspace(0, 3, 101)
yy = exacte(tt)

plt.plot(tt,yy);

Le schéma de Crank-Nicolson construit la suite \((u_n)_{n\in \mathbb{N}}\) selon la relation \( u_{n+1}=u_n-5h(u_{n+1}+u_n), \qquad \forall n \in \mathbb{N}, \) pour \(h>0\) fixé. On obtient une formule de récurrence rendue explicite par un calcul élémentaire: \( u_{n+1}=-5hu_{n+1}-5hu_n+u_n \) d’où \( u_{n+1}=\frac{1-5h}{1+5h}u_n. \) Par récurrence on obtient \( u_{n}=\left(\frac{1-5h}{1+5h}\right)^nu_0. \) Il s’agit d’une suite géométrique de raison \( r=\frac{1-5h}{1+5h}. \)

Pour tout \(h>0\) on a \( r=\frac{1-5h}{1+5h}=1-\frac{10h}{1+5h} \) et \(-1<1-\frac{10h}{1+5h}<1. \) Ce schéma est donc inconditionnellement A-stable car \(|u_{n+1}|=|r^{n+1}u_0|\le|u_0|\).

Le schéma génère une suite positive ssi \( 1-\frac{10h}{1+5h}>0 \) c’est-à-dire ssi \( h<\frac{1}{5}. \)

import matplotlib.pyplot as plt
import numpy as np

# Réglage de la taille de la police
plt.rcParams.update({'font.size': 12})

N = 6
tt = np.linspace(0, 1, N+1)
y0 = 1

exacte = lambda t: y0 * np.exp(-10*t)

phi = lambda t, y: -y**5
def schema(phi, tt, y0):
    h = tt[1] - tt[0]
    uu = np.empty_like(tt)
    uu[0] = y0
    for i in range(len(tt)-1):
        uu[i+1] = (1 - 5*h) / (1 + 5*h) * uu[i]
    return uu


yy = exacte(tt)
uu = schema(phi, tt, y0)

plt.figure(figsize=(8, 6))
plt.plot(tt, yy, label="Exacte")
plt.plot(tt, uu, "o", label="Schéma")
plt.title("Comparaison de la solution exacte et du schéma")
plt.grid()
plt.legend();

Exercice : A-stabilité d’un schéma

Soit le problème de Cauchy:

\(\begin{cases} y'(t)+\dfrac{\sqrt{y(t)}}{2}=0, & \forall t \in \mathbb{R}^+,\\ y(0)=y_0>0. \end{cases}\)

Pour \(h>0\) fixé, soit le schéma numérique défini par la suite \((u_n)_{n\in \mathbb{N}}\) suivante

\(\frac{u_{n+1}-u_n}{h}+\frac{u_{n+1}}{2\sqrt{u_n}}=0, \qquad \forall n \in \mathbb{N},\)

1. Expliciter l’expression de \(u_{n+1}\) en fonction de \(u_n\).
2. Montrer que la suite \((u_n)_{n\in \mathbb{N}}\) est une suite positive, décroissante et convergente vers \(0\).

Correction

C’est un problème de Cauchy du type \(\begin{cases} y'(t)=\varphi(t,y(t)), & \forall t \in \mathbb{R},\\ y(0)=y_0>0, \end{cases}\) avec \(\varphi(t,y)=g(y)=-\frac{\sqrt{y}}{2}\).

Pour \(h>0\) fixé on obtient une formule de récurrence rendue explicite par un calcul élémentaire: \( u_{n+1}=\frac{u_n}{1 + \frac{h}{2\sqrt{u_n}}}=\frac{2(u_n)^{3/2}}{2\sqrt{u_n}+h}, \qquad \forall n \in \mathbb{N}. \)

On étudie la suite \(\begin{cases} u_0>0,\\ u_{n+1}=\frac{2(u_n)^{3/2}}{2\sqrt{u_n}+h}, \qquad \forall n \in \mathbb{N}, \end{cases}\) pour \(h>0\) fixé.

Par récurrence on montre que si \(u_0>0\) alors \(u_n>0\) pour tout \(n\) \(\in\) \(\mathbb{N}\). De plus, \(\frac{u_{n+1}}{u_{n}}=\frac{1}{1 + \frac{h}{2\sqrt{u_n}}}<1\) pour tout \(n\) \(\in\) \(\mathbb{N}\): la suite est monotone décroissante. On a alors une suite décroissante et bornée par zéro, donc elle converge. Soit \(\ell\) la limite de cette suite, alors \(\ell\ge0\) et \(\ell= \frac{2\ell^{3/2}}{2\sqrt{\ell}+h}\) donc \(\ell=0\).

    import matplotlib.pyplot as plt
import numpy as np

# Réglage de la taille de la police
plt.rcParams.update({'font.size': 16})

t0 = 0
y0 = 1
N = 10
tt = np.linspace(t0, 10, N+1)


phi = lambda t, y: -0.5 * np.sqrt(y)
exacte = lambda t: np.where(t >= 4 * np.sqrt(y0), 0, (np.sqrt(y0) - t / 4) ** 2)


def schema(phi, tt, y0):
    h = tt[1] - tt[0]
    uu = [y0]
    for i in range(len(tt)-1):
        uu.append(2 * (uu[i]) ** (3/2) / (h + 2 * np.sqrt(uu[i])))
    return uu

yy = exacte(tt)
uu = schema(phi, tt, y0)

plt.figure(figsize=(8, 6))
plt.plot(tt, yy, label="Exacte")
plt.plot(tt, uu, "o", label="Schéma")
plt.title("Comparaison de la solution exacte et du schéma")
plt.grid()
plt.legend();

Exercice : A-stabilité d’un schéma

Soit le problème de Cauchy:

\(\begin{cases} y'(t)+y^5(t)=0, & \forall t \in \mathbb{R}^+,\\ y(0)=y_0>0. \end{cases}\)

1. Montrer qu’il existe une unique solution globale \(y\) \(\in\) \(\mathscr{C}^\infty(\mathbb{R}^+,\mathbb{R}^+)\).
2. Soit le schéma numérique défini par la suite \((u_n)_{n\in \mathbb{N}}\) suivante \( \frac{u_{n+1}-u_n}{h}+u_{n+1}u_n^4=0, \qquad \forall n \in \mathbb{N}, \) pour \(h>0\) fixé. Expliciter l’expression de \(u_{n+1}\) en fonction de \(u_n\).
3. Montrer que la suite \((u_n)_{n\in \mathbb{N}}\) est une suite décroissante et convergente vers \(0\) pour tout \(h>0\) fixé.

Correction

C’est un problème deCauchy du type \(\begin{cases} y'(t)=\varphi(t,y(t)), & \forall t \in \mathbb{R},\\ y(0)=y_0>0, \end{cases}\) avec \(\varphi(t,y)=g(y)=-y^5\).

On montre d’abord qu’il existe une et une seule solution locale (ie sur \([0;T]\)) de classe \(\mathscr{C}^1([0,T],\mathbb{R})\).
On montre ensuite que cette solution est de classe \(\mathscr{C}^\infty([0,T],\mathbb{R})\).
On montre enfin que la solution admet un prolongement sur \(\mathbb{R}\). + Comme \(g\) \(\in\) \(\mathscr{C}^1(\mathbb{R}^+,\mathbb{R}^+)\), d’après le théorème de Cauchy-Lipschitz il existe \(T>0\) et une unique solution \(y\) \(\in\) \(\mathscr{C}^1([0,T],\mathbb{R})\). + Par récurrence, en exploitant l’EDO et la régularité de \(g\), on grimpe en régularité sur \(y\) ainsi \(y\) \(\in\) \(\mathscr{C}^\infty([0,T],\mathbb{R})\). + La fonctionne nulle est solution de l’EDO (mais non du problème de Cauchy donné). Comme \(y_0>0\), par l’unicité de la solution du problème de Cauchy on a \(y(t)>0\) pour tout \(t\) \(\in\) \([0,T]\) (car deux trajectoires ne peuvent pas se croiser). De plus, \(y\) est décroissante, ainsi la solution est bornée (\(y(t)\) \(\in\) \(]0,y_0[\)). On en déduit par le théorème des extrémités que la solution \(y\) admet un prolongement sur \(\mathbb{R}^+\) solution de l’EDO.

Pour en calculer la solution, on remarque qu’il s’agit d’une EDO à variables séparables. L’unique solution du problème de Cauchy est \( y(t)=y_0(4ty_0^4+1)^{-1/4} \)

import matplotlib.pyplot as plt
import numpy as np
plt.rcParams.update({'font.size': 12})

t0 = 0
y0 = 1
exacte = lambda t: y0 * (4 * t * y0 ** 4 + 1) ** (-0.25)

tt = np.linspace(t0, 10, 101)
yy = exacte(tt)

plt.plot(tt, yy)
plt.grid()

Pour \(h>0\) fixé on obtient une formule de récurrence rendue explicite par un calcul élémentaire: \( u_{n+1}=\frac{u_n}{1 + u_n^4 h}, \qquad \forall n \in \mathbb{N}. \)

On étudie la suite \(\begin{cases} u_0=y_0>0,\\ u_{n+1}=\frac{u_n}{1 + u_n^4 h}, \qquad \forall n \in \mathbb{N}, \end{cases}\) pour \(h>0\) fixé.

Par récurrence on montre que si \(u_0>0\) alors \(u_n>0\) pour tout \(n\) \(\in\) \(\mathbb{N}\). De plus, \(\frac{u_{n+1}}{u_{n}}=\frac{1}{1 + u_n^4 h}<1\) pour tout \(n\) \(\in\) \(\mathbb{N}\): la suite est monotone décroissante. On a alors une suite décroissante et bornée par zéro, donc elle converge. Soit \(\ell\) la limite de cette suite, alors \(\ell\ge0\) et \(\ell= \frac{\ell}{1 + \ell^4 h}\) donc \(\ell=0\). Ce schéma est donc inconditionnellement A-stable.

    import matplotlib.pyplot as plt
import numpy as np

# Réglage de la taille de la police
plt.rcParams.update({'font.size': 16})

t0 = 0
y0 = 1
N = 50
tt = np.linspace(t0, 10, N+1)


phi = lambda t, y: -y ** 5
exacte = lambda t: y0 * (4 * t * y0 ** 4 + 1) ** (-0.25)


def schema(phi, tt, y0):
    h = tt[1] - tt[0]
    uu = np.empty_like(tt)
    uu[0] = y0
    for i in range(len(tt) - 1):
        uu[i + 1] = uu[i] / (1 + h * uu[i] ** 4)
    return uu

yy = exacte(tt)
uu = schema(phi, tt, y0)

plt.figure(figsize=(8, 6))
plt.plot(tt, yy, label="Exacte")
plt.plot(tt, uu, "o", label="Schéma")
plt.title("Comparaison de la solution exacte et du schéma")
plt.grid()
plt.legend();

Exercice : schéma d’Euler explicite et stabilité

Soit le problème de Cauchy:

\(\begin{cases} y'(t)+\sin(y(t))=0, & \forall t \in \mathbb{R},\\ y(0)=y_0>0. \end{cases}\)

1. Montrer qu’il existe une unique solution globale \(y\) \(\in\) \(\mathscr{C}^\infty(\mathbb{R},\mathbb{R})\).
2. Écrire le schéma le schéma d’Euler explicite pour ce problème.
3. Montrer que la suite \((u_n)_{n\in \mathbb{N}}\) construite par ce schéma vérifie \( |u_{n+1}|\le |u_n|+h, \qquad \forall n \in \mathbb{N},\) où \(h>0\) est le pas de la suite.
4. En déduire que

\( |u_{n}|\le |u_0|+nh, \qquad \forall n \in \mathbb{N}. \)

Correction

C’est un problème deCauchy du type \(\begin{cases} y'(t)=\varphi(t,y(t)), & \forall t \in \mathbb{R},\\ y(0)=y_0>0, \end{cases}\) avec \(\varphi(t,y)=g(y)=-\sin(y)\).

On montre d’abord qu’il existe une et une seule solution locale (ie sur \([0;T]\)) de classe \(\mathscr{C}^1([0,T],\mathbb{R})\).
On montre ensuite que cette solution est de classe \(\mathscr{C}^\infty([0,T],\mathbb{R})\).
On montre enfin que la solution admet un prolongement sur \(\mathbb{R}\). + Comme \(g\) \(\in\) \(\mathscr{C}^1(\mathbb{R},\mathbb{R})\), d’après le théorème de Cauchy-Lipschitz il existe \(T>0\) et une unique solution \(y\) \(\in\) \(\mathscr{C}^1([-T,T],\mathbb{R})\). + Par récurrence, en exploitant l’EDO et la régularité de \(g\), on grimpe en régularité sur \(y\) ainsi \(y\) \(\in\) \(\mathscr{C}^\infty([-T,T],\mathbb{R})\). + Toutes les fonctions constante \(y(t)=k\pi\) pour \(k\in\mathbb{Z}\) sont solutions de l’équation différentielle car \(g(k\pi)=0\). Pour \(y_0\) donné, soit \(\tilde k\in\mathbb{Z}\) tel que \(y_0\in[\tilde k\pi;(\tilde k+1)\pi]\); par l’unicité de la solution du problème de Cauchy on a \(y(t)\in[\tilde k\pi;(\tilde k+1)\pi]\) pour tout \(t\) \(\in\) \([-T,T]\) (car deux trajectoires ne peuvent pas se croiser), i.e. la solution est bornée. On en déduit par le théorème des extrémités que la solution \(y\) admet un prolongement sur \(\mathbb{R}\) solution de l’EDO.

Soit \(h>0\) fixé et \(t_n=nh\) pour tout \(n\in\mathbb{Z}\). Le schéma d’Euler explicite pour l’EDO donnée construit la suite \( u_{n+1}=u_n - h \sin(u_n), \qquad \forall n \in \mathbb{N}. \)

Comme \(|a+b|\le |a|+|b|\) et comme \(|-\sin(x)|\le1\) pour tout \(x\in\mathbb{R}\), on conclut que \( |u_{n+1}|=|u_n - h \sin(u_n)| \le |u_n| + |h \sin(u_n)| \le |u_n| + h \) pour \(h>0\) fixé.

Par récurrence: \(|u_{n+1}|\le |u_n| + h \le |u_{n-1}| + 2h \le\dots\le |u_0| + (n+1)h\).

t0 = 0 
y0 = 1

# SOLUTION EXACTE POUR y(0) = 1

import sympy as sym
sym.init_printing()

t  = sym.Symbol('t')
y  = sym.Function('y')

right = -sym.sin(y(t)) # phi(t,y(t)) avec sympy
edo = sym.Eq( sym.diff(y(t),t) , right )
display(edo)

solgenLIST = sym.dsolve(edo,y(t))
display(solgenLIST)

solgen = solgenLIST[1]
display(solgen)

equation = sym.Eq( y0, solgen.rhs.subs(t,t0))
display(equation)

consts = sym.solve( equation , dict=True)[0]
display(consts)

solpar = solgen.subs(consts)
display(solpar)

func = sym.lambdify(t,solpar.rhs,'numpy')

\(\displaystyle \frac{d}{d t} y{\left(t \right)} = - \sin{\left(y{\left(t \right)} \right)}\)

\(\displaystyle \left[ y{\left(t \right)} = - \operatorname{acos}{\left(\frac{- C_{1} - e^{2 t}}{C_{1} - e^{2 t}} \right)} + 2 \pi, \ y{\left(t \right)} = \operatorname{acos}{\left(\frac{- C_{1} - e^{2 t}}{C_{1} - e^{2 t}} \right)}\right]\)

\(\displaystyle y{\left(t \right)} = \operatorname{acos}{\left(\frac{- C_{1} - e^{2 t}}{C_{1} - e^{2 t}} \right)}\)

\(\displaystyle 1 = \operatorname{acos}{\left(\frac{- C_{1} - 1}{C_{1} - 1} \right)}\)

\(\displaystyle \left\{ C_{1} : \frac{-1 + \cos{\left(1 \right)}}{\cos{\left(1 \right)} + 1}\right\}\)

\(\displaystyle y{\left(t \right)} = \operatorname{acos}{\left(\frac{- e^{2 t} - \frac{-1 + \cos{\left(1 \right)}}{\cos{\left(1 \right)} + 1}}{- e^{2 t} + \frac{-1 + \cos{\left(1 \right)}}{\cos{\left(1 \right)} + 1}} \right)}\)

# SOLUTION APPROCHEE PAR EULER EXPLICITE POUR y(0) = 1

import numpy as np

phi = lambda t, y: -np.sin(y)

def schema(phi, tt, y0):
    h = tt[1] - tt[0]
    uu = np.empty_like(tt)
    uu[0] = y0
    for i in range(len(tt) - 1):
        uu[i + 1] = uu[i]+h*phi(tt[i],uu[i])
    return uu

# COMPARAISONS
N = 100
tt = np.linspace(t0,10,N+1)

yy = func(tt)
uu = schema(phi, tt, y0)

import matplotlib.pyplot as plt 
plt.rcParams.update({'font.size': 16})

plt.figure(figsize=(10,7))
plt.plot(tt, yy, label="Exacte")
plt.plot(tt, uu, 'o', label="Schéma")
plt.legend();

Exercice : A-stabilité d’un schéma

Un modèle pour la diffusion d’une épidémie se base sur l’hypothèse que sa vitesse de propagation est proportionnelle au nombre d’individus infectés et au nombre d’individus sains.

Si on note \(I(t)\ge0\) le nombre d’individus infectés à l’instant \(t\ge0\) et \(A>0\) le nombre d’individus total, il existe une constante \(k\) \(\in\) \(\mathbb{R}^+\) telle que \(I'(t)=k I(t)(A-I(t))\).

1. Montrer qu’il existe \(T>0\) et une unique solution \(I\) \(\in\) \(\mathscr{C}^\infty([0,T])\) au problème de Cauchy: \( \begin{cases} I'(t)=k I(t)(A-I(t)),\\ I(0)=I_0>0. \end{cases} \)
   
2. Montrer que si \(0<I_0<A\) alors \(0<I(t)<A\) pour tout \(t>0\).
   
3. Montrer que si \(0<I_0<A\) alors \(I(t)\) est croissante sur \(\mathbb{R}^+\).
4. Soit \(0<I_0<A\). On considère le schéma semi-implicite \(\frac{I_{n+1}-I_n}{h}=kI_n(A-I_{n+1}).\) Montrer que \(I_n\to A\) lorsque \(n\to+\infty\) indépendamment du pas \(h>0\) fixé.
5. Calculer la solution exacte

Correction

C’est un problème deCauchy du type \(\begin{cases} I'(t)=\varphi(t,I(t)), & \forall t \in \mathbb{R}^+,\\ I(0)=I_0>0, \end{cases}\) avec \(\varphi(t,I(t))=g(I(t))=kI(t)(A-I(t))\).

• Comme \(g\) \(\in\) \(\mathscr{C}^1(\mathbb{R},\mathbb{R})\), d’après le théorème de Cauchy-Lipschitz il existe \(T>0\) et une unique solution \(I\) \(\in\) \(\mathscr{C}^1([0,T],\mathbb{R})\).
• Par récurrence, en exploitant l’EDO et la régularité de \(g\), on grimpe en régularité sur \(I\) ainsi \(I\) \(\in\) \(\mathscr{C}^\infty([0,T],\mathbb{R})\).
• Puisque la fonction nulle et la fonction constante \(I(t)=A\) sont solutions de l’équation différentielle, si \(0<I_0<A\) alors \(0<I(t)<A\) pour tout \(t\) \(\in\) \([0,T]\) (car, par l’unicité de la solution du problème de Cauchy, deux trajectoires ne peuvent pas se croiser).
• Puisque \(I'(t)=kI(t)(A-I(t))\), si \(0<I_0<A\) alors \(I\) est croissante pour tout \(t\) \(\in\) \([0,T]\). On en déduit par le théorème des extrémités que la solution \(I\) admet un prolongement sur \(\mathbb{R}^+\) solution de l’EDO et que \(I\) est croissante pour tout \(t\) \(\in\) \(\mathbb{R}^+\).

Soit \(0<I_0<A\). On considère le schéma semi-implicite \( \frac{I_{n+1}-I_n}{h}=kI_n(A-I_{n+1}) \) pour \(h>0\) fixé. On obtient une formule de récurrence rendue explicite par un calcul élémentaire: \( I_{n+1}=\frac{1+kAh}{1+kI_nh}I_n. \) Si \(0<I_0<A\) alors + \(I_n>0\) quelque soit \(n\); + \(I_n\) est majorée par \(A\) car \( I_{n+1} \le A \quad\iff\quad (1 + k A h ) I_n \le (1 + k I_n h) A \quad\iff\quad I_n \le A \) donc par récurrence \(I_{n+1} \le A\) quelque soit \(n\); + si \(I_n\to \ell\) alors \(\ell=\frac{1+k A h}{1+k\ell h}\ell\) donc \(\ell=0\) ou \(\ell=A\); + \(I_n\) est une suite monotone croissante (encore par récurrence on montre que \(|I_{n+1}|\ge|I_n|\ge\dots\ge|I_0|\));

donc \(I_n\to A\) lorsque \(n\to+\infty\) indépendamment du pas \(h>0\) choisi.

On a déjà observé qu’il y a deux solutions constantes de l’EDO: la fonction \(I(t)\equiv0\) et la fonction \(I(t)\equiv A\).

Pour chercher toutes les solutions non constantes on remarque qu’il s’agit d’une EDO à variables séparables donc on a \(\begin{aligned} % I'(t)=k I(t)(5000-I(t))\\ % \frac{I'(t)}{I(t)(5000-I(t))}&=k\\ % \frac{\mathrm{d}\,I}{I(5000-I)}&=k\mathrm{d}\,t\\ % \int\frac{1}{I(5000-I)}\mathrm{d}\,I&=k\int\mathrm{d}\,t\\ % \int\frac{1}{I}\mathrm{d}\,I-\int\frac{1}{5000-I}\mathrm{d}\,I&=5000k\int\mathrm{d}\,t\\ % \ln(I)+\ln(5000-I)&=5000kt+c\\ % \ln\frac{I}{5000-I}&=5000kt+c\\ % \frac{I}{5000-I}&=De^{5000kt}\\ % I(t)&=\frac{5000De^{5000 k t}}{1+De^{5000 k t}}\\ I(t)&=\frac{A}{De^{-A k t}+1} \end{aligned}\)

La valeur numérique de la constante d’intégration \(D\) est obtenue grâce à la CI: \( D=\frac{A-I_0}{I_0} \)

Exemple avec \(A=5000\), \(I_0=160\), \(k=\frac{\ln(363/38)}{35000}\) et \(h=1\):

import numpy as np
import matplotlib.pyplot as plt

# DONNEES
A = 5000
I0 = 160
k = np.log(363/38) / 35000

D = (A - I0) / I0
phi = lambda t, I: k * I * (A - I) # inutile pour notre schema


# SOLUTION EXACTE
exacte = lambda t: A / (D * np.exp(-A * k * t) + 1)

# SCHEMA
def schema(phi, tt, y0):
    uu = np.empty_like(tt)
    uu[0] = y0
    for i in range(len(tt) - 1):
        uu[i + 1] = (1 + k * A * h) / (1 + k * uu[i] * h) * uu[i]
    return uu


# MAIN
h = 1
tt = np.arange(0, 30, h)

yy = exacte(tt)
uu = schema(phi, tt, I0)

plt.plot(tt, yy, label="Exacte")
plt.plot(tt, uu, 'o', label="Schéma")
plt.legend();

Exercice : A-stabilité du \(\vartheta\)-schéma

Pour calculer une approximation de la solution du problème de Cauchy

\(\begin{cases} y(t_0)=y_0,\\ y'(t)=\varphi(t,y(t)). \end{cases}\)

considérons la \(\vartheta\)-méthode définie comme suit:

\(\begin{cases} u_0=y_0,\\ u_{n+1}=u_{n}+h\Big(\vartheta\varphi(t_{n+1},u_{n+1})+(1-\vartheta)\varphi(t_{n},u_{n})\Big). \end{cases}\)

Ici \(\vartheta\) est un paramètre dans l’intervalle \([0;1]\).

L’objectif de cet exercice est d’étudier la A-stabilité de cette méthode en fonction du paramètre \(\vartheta\). Pour cela considérons l’EDO obtenue avec \(\varphi(t,y(t))=-\beta y(t)\) avec \(\beta>0\).

Donner une condition sur \(h\) pour que la \(\vartheta\)-méthode soit A-stable. Pour quelles valeurs de \(\vartheta\) la méthode est inconditionnellement A-stable?

Remarque : pour \(\vartheta=0\) on retrouve la méthode d’Euler explicite, pour \(\vartheta=1\) la méthode d’Euler implicite, pour \(\vartheta=\frac{1}{2}\) la méthode de Crank-Nicholson. On voit donc que pour certaines valeurs de \(\vartheta\) la méthode est inconditionnellement A-stable et pour d’autres elle ne l’est pas.

Correction

\( u_{n+1}=u_n-\beta h(\vartheta u_{n+1}+(1-\vartheta)u_n) \iff u_{n+1} =\frac{1-\beta h (1-\vartheta)}{1+\beta h \vartheta} u_n =\frac{(1+\beta h\vartheta)-\beta h }{1+\beta h \vartheta} u_n \) Par induction \( u_n=\left(1-\frac{\beta h }{1+\beta h \vartheta}\right)^n y_0 \) Il s’agit d’une suite géométrique de raison \(q=1-\frac{\beta h }{1+\beta h \vartheta}\).

La suite tends à zéro ssi \(|q|<1\) (la convergence est monotone ssi \(0\le q<1\)): \(-1<1-\frac{\beta h }{1+\beta h \vartheta}<1 \iff 0<\frac{\beta h }{1+\beta h \vartheta}<2 \iff (1-2\vartheta)\beta h < 2 \) Pour \(\vartheta\ge \frac{1}{2}\) la méthode est inconditionnellement A-stable,
pour \(\vartheta< \frac{1}{2}\) la méthode est A-stable ssi \(h<\frac{2}{(1-2\vartheta)\beta}\).

Exercice : Déduction d’un schéma Predictor-Corrector

Écrire le schéma predictor-corrector basé sur AM-3 AB-3. Attention à bien initialiser la suite récurrente.

Correction

• AB-1 \(u_{n+1}=u_n+h\varphi(t_n,u_n)\)
• AB-2 \(u_{n+1}=u_n+\frac{h}{2}\left(3\varphi(t_n,u_n)-\varphi(t_{n-1},u_{n-1})\right)\)
• AB-3 \(u_{n+1}=u_n+\frac{h}{12}\left(23\varphi(t_n,u_n)-16\varphi(t_{n-1},u_{n-1})+5\varphi(t_{n-2},u_{n-2})\right)\)
• AM-3 \(u_{n+1}=u_n+\frac{h}{24}\left(9\varphi(t_{n+1},u_{n+1})+19\varphi(t_n,u_n)-5\varphi(t_{n-1},u_{n-1})+\varphi(t_{n-2},u_{n-2})\right)\) donc \(\begin{cases} u_0=y_0,\\ u_1=u_0+h\varphi(t_0,u_0)\\ u_{2}=u_1+\frac{h}{2}\left(3\varphi(t_1,u_1)-\varphi(t_0,u_0)\right) \\ \\ \tilde u=u_n+\frac{h}{12}\left(23\varphi(t_n,u_n)-16\varphi(t_{n-1},u_{n-1})+5\varphi(t_{n-2},u_{n-2})\right)\\ u_{n+1}=u_n+\frac{h}{24}\left(9\varphi(t_{n+1},\tilde u)+19\varphi(t_n,u_n)-5\varphi(t_{n-1},u_{n-1})+\varphi(t_{n-2},u_{n-2}\right) \end{cases} \) AM-3 est un schéma implicite à 3 pas d’ordre 4,
  AB-3 est un schéma explicite à 3 pas d’ordre 3,
  le schéma predictor-corrector est donc un schéma à 3 pas d’ordre 4.

Exercice : Ordre d’un schéma

Une méthode numérique à un pas a été utilisée pour résoudre une équation différentielle avec condition initiale. Les résultats obtenus par cette méthode en prenant des pas de temps \(h = 0.1\), \(h = 0.05\) et \(h = 0.025\) sont donnés dans le tableau suivant (remarque: une valeur sur deux est affichée pour \(h = 0.05\) et une valeur sur quatre est affichée pour \(h = 0.025\))

\(\begin{array}{|c|c|c|c|} \hline t_i &y_i\text{ pour }h=0.1 &y_i\text{ pour }h=0.05 &y_i\text{ pour }h=0.025\\ \hline 1.0 &0.500000 &0.500000 &0.500000\\ 1.1 &0.512084 &0.512242 &0.512280\\ 1.2 &0.511698 &0.512101 &0.512196\\ 1.3 &0.500927 &0.501559 &0.501704\\ 1.4 &0.482686 &0.483447 &0.483619\\ 1.5 &0.459861 &0.460633 &0.460804\\ \hline \end{array}\)

En calculant le rapport des erreurs et sachant que \(y(1.5)=0.460857\), déterminer l’ordre de la méthode numérique utilisée.

Correction

La méthode utilisée est d’ordre \(2\), en effet, pour \(t=1.5\) nous avons \( \frac{E(h = 0. 05)}{E(h = 0. 025)}=\frac{|0.460633 - 0.460857|}{|0.460804 - 0.460857|}= 4.226\simeq2^2. \) Si on a accès à polyfit on peut tracer la courbe d’erreur en echèlle logaritmique et estimer la pente:

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 12})


H = [0.1, 0.05, 0.025]
# Calcul des erreurs pour chaque valeur de H au point final
err = []
y = 0.460857
err.append(abs(y - 0.459861))
err.append(abs(y - 0.460633))
err.append(abs(y - 0.460804))

# Calcul de l'ordre de convergence
order = np.polyfit(np.log(H), np.log(err), 1)[0]
print('Ordre de convergence :', round(order, 2))

plt.figure(figsize=(12, 6))
# Tracé du graphique en double-log pour visualiser la convergence
plt.subplot(1, 2, 1)
plt.loglog(H, err, 'r-o')
plt.grid()

# Tracé du graphique des logarithmes de l'erreur en fonction du logarithme du pas de temps
plt.subplot(1, 2, 2)
plt.plot(np.log(H), np.log(err), 'r-o');

Ordre de convergence : 2.12

Exercice : Interpolation, Quadrature et EDO

• Soit \(f\) une fonction de classe \(\mathcal{C}^1([-1,1])\). Écrire le polynôme \(p\in\mathbb{R}_2[\tau]\) qui interpole \(f\) aux points \(-1\), \(0\) et \(1\).

• Construire une méthode de quadrature comme suit: \( \int_{0}^{1}f(\tau)\mathrm{d}\tau \approx \int_{0}^{1}p(\tau)\mathrm{d}\tau. \) NB: on intègre sur \([0,1]\) mais on interpole en \(-1\), \(0\) et \(1\).

• À l’aide d’un changement de variable affine entre l’intervalle \([0,1]\) et l’intervalle \([a,b]\), en déduire une formule de quadrature pour l’intégrale \( \int_{a}^{b}f(x) \mathrm{d}x \) lorsque \(f\) est une fonction de classe \(\mathcal{C}^1([2a-b,b])\).
  Remarque: \([2a-b,b]=[a-(b-a),a+(b-a)]\)

• Considérons le problème de Cauchy: trouver $ y $ tel que \( \begin{cases} y'(t) = \varphi(t,y(t)), &\forall t \in [t_0,T],\\ y(t_0) = y_0, \end{cases} \) dont on suppose l’existence d’une unique solution \(y\).

  • On subdivise l’intervalle \([t_0;T]\) en \(N\) intervalles \([t_{n};t_{n+1}]\) de largeur \(h=\dfrac{T-t_0}{N}\) avec \(t_n=t_0+nh\) pour \(n=0,\dots,N\).
  • Utiliser la formule obtenue au point 3 pour approcher l’intégrale \( \int_{t_n}^{t_{n+1}}\varphi(t,y(t))\mathrm{d}t. \)
  • En déduire un schéma à deux pas implicite pour l’approximation de la solution du problème de Cauchy.

Correction

On cherche les coefficients \(\alpha\), \(\beta\) et \(\gamma\) du polynôme \(p(\tau)=\alpha+\beta \tau+\gamma \tau^2\) tels que \( \begin{cases} p(-1)=f(-1),\\ p(0) =f(0),\\ p(1) =f(1), \end{cases} \qquad\text{c'est à dire}\qquad \begin{cases} \alpha-\beta+\gamma=f(-1),\\ \alpha=f(0),\\ \alpha+\beta+\gamma=f(1). \end{cases} \) Donc \(\alpha=f(0)\), \(\beta=\frac{f(1)-f(-1)}{2}\) et \(\gamma=\frac{f(1)-2f(0)+f(-1)}{2}\).

import sympy as symb
symb.init_printing()

f_0 = symb.Symbol('f(0)')
f_1 = symb.Symbol('f(1)')
f_m1 = symb.Symbol('f(-1)')
t = symb.Symbol('t')

p = symb.interpolate([(1,f_1),(0,f_0),(-1,f_m1)], t).factor(t)
display(symb.Eq(symb.Symbol('p'), p))
q = (symb.integrate(p,(t,0,1)).simplify())
display(symb.Eq( symb.Integral(symb.Symbol('p(t)'), (t, 0, 1)) , q ))

\(\displaystyle p = \frac{2 f(0) + t^{2} \left(f(-1) - 2 f(0) + f(1)\right) + t \left(- f(-1) + f(1)\right)}{2}\)

\(\displaystyle \int\limits_{0}^{1} p(t)\, dt = - \frac{f(-1)}{12} + \frac{2 f(0)}{3} + \frac{5 f(1)}{12}\)

On en déduit la méthode de quadrature \( \int_{0}^{1}f(\tau)\mathrm{d}\tau\approx\int_{0}^{1}p(\tau)\mathrm{d}\tau=\alpha+\frac{\beta}{2}+\frac{\gamma}{3}=\frac{-f(-1)+8f(0)+5f(1)}{12}. \)

Soit \(x=m\tau+q\), alors \( \int_{a}^{b} f(x)\mathrm{d}x=m\int_{0}^{1}f(m\tau+q)\mathrm{d}\tau \quad\text{avec}\quad \begin{cases} a=q,\\ b=m+q, \end{cases} \quad\text{i.e.}\quad \begin{cases} m=b-a,\\ q=a, \end{cases} \) d’où le changement de variable \(x=(b-a)\tau+a\). On en déduit la formule de quadrature \( \int_{a}^{b}\!\!\!\!f(x)\mathrm{d}x=(b-a)\int_{0}^{1}f\left((b-a)\tau+a\right)\mathrm{d}\tau\approx(b-a)\frac{-f(2a-b)+8f(a)+5f(b)}{12}. \)

On pose \(a=t_n\) et \(b=t_{n+1}\) d’où la formule de quadrature \(\int_{t_{n}}^{t_{n+1}}\!\!\!\!f(t)\mathrm{d}t\approx(t_{n+1}-t_{n}) \frac{-f(2t_{n}-t_{n+1})+8f(t_{n})+5f(t_{n+1})}{12}=h \frac{-f(t_{n-1})+8f(t_{n})+5f(t_{n+1})}{12}.\)

En utilisant la formule de quadrature pour l’intégration de l’EDO \(y'(t)=\varphi(t,y(t))\) entre \(t_n\) et \(t_{n+1}\) on obtient \( y(t_{n+1})=y(t_n)+\int_{t_n}^{t_{n+1}} \varphi(t,y(t))\mathrm{d}t\approx h \dfrac{-\varphi(t_{n-1},y(t_{n-1}))+8\varphi(t_{n},y(t_{n}))+5\varphi(t_{n+1},y(t_{n+1}))}{12}. \)

Si on note \(u_n\) une approximation de la solution \(y\) au temps \(t_n\), on obtient le schéma à deux pas implicite suivant: \( \begin{cases} u_0=y(t_0)=y_0,\\ u_1\text{ à définir},\\ u_{n+1}=u_n+ h \dfrac{-\varphi(t_{n-1},u_{n-1})+8\varphi(t_{n},u_{n})+5\varphi(t_{n+1},u_{n+1})}{12} & n=1,2,\dots N-1 \end{cases} \) On peut utiliser une pédiction d’Euler explicite pour initialiser \(u_1\): \( \begin{cases} u_0=y(t_0)=y_0,\\ u_1=u_0+h\varphi(t_0,u_0),\\ u_{n+1}=u_n+ h \dfrac{-\varphi(t_{n-1},u_{n-1})+8\varphi(t_{n},u_{n})+5\varphi(t_{n+1},u_{n+1})}{12} & n=1,2,\dots N-1 \end{cases} \)

Exercice : Interpolation, Quadrature et EDO

1. Soit \(h>0\) et \(f\colon [a-h,a+h] \to\mathbb{R}\) une fonction de classe \(\mathscr{C}^1([a-h,a+h])\). Écrire le polynôme \(p\in\mathbb{R}_1[x]\) qui interpole \(f\) aux points \(a-h\) et \(a\), i.e. l’équation de la droite \(p\in\mathbb{R}_1[x]\) qui passe par les deux points \((a-h,f(a-h))\) et \((a,f(a))\).
2. Construire une méthode de quadrature comme suit : \( \int_{a}^{a+h}f(x)\mathrm{d}x \approx \int_{a}^{a+h}p(x)\mathrm{d}x. \) NB: on intègre sur \([a,a+h]\) mais on interpole en \(a-h\) et \(a\).
3. Considérons le problème de Cauchy: trouver \(y \colon [t_0,T]\subset \mathbb{R} \to \mathbb{R}\) tel que \( \begin{cases} y'(t) = \varphi(t,y(t)), &\forall t \in [t_0,T],\\ y(t_0) = y_0, \end{cases} \) dont on suppose l’existence d’une unique solution \(y\). On subdivise l’intervalle \([t_0;T]\) en \(N\) intervalles \([t_{n};t_{n+1}]\) de largeur \(h=\dfrac{T-t_0}{N}\) avec \(t_n=t_0+nh\) pour \(n=0,\dots,N\).
   Utiliser la formule obtenue au point 2 pour approcher l’intégrale \( \int_{t_n}^{t_{n+1}}\varphi(t,y(t))\mathrm{d}t. \) En déduire un schéma à deux pas explicite pour l’approximation de la solution du problème de Cauchy.

Correction

\(p(x)=\dfrac{f(a)-f(a-h)}{a-(a-h)}(x-a)+f(a)=\dfrac{f(a)-f(a-h)}{h}(x-a)+f(a)\).

import sympy as symb
symb.init_printing()

a = symb.Symbol('a')
h = symb.Symbol('h')
f_a = symb.Symbol('f(a)')
f_amh = symb.Symbol('f(a-h)')

p = symb.interpolate([(a,f_a),(a-h,f_amh)], t).factor(t)
display(symb.Eq(symb.Symbol('p'), p))
q = (symb.integrate(p,(t,a,a+h)).simplify())
display(symb.Eq( symb.Integral(symb.Symbol('p(t)'), (t, a, a+h)) , q ))

\(\displaystyle p = \frac{- a f(a) + a f(a-h) + f(a) h + t \left(f(a) - f(a-h)\right)}{h}\)

\(\displaystyle \int\limits_{a}^{a + h} p(t)\, dt = \frac{h \left(3 f(a) - f(a-h)\right)}{2}\)

On en déduit la méthode de quadrature

\(\begin{aligned} \int_{a}^{a+h}f(x)\mathrm{d}x &\approx \int_{a}^{a+h}p(x)\mathrm{d}x \\&= \dfrac{f(a)-f(a-h)}{h}\left[ \frac{(x-a)^2}{2} \right]_a^{a+h}+f(a)\left[ x \right]_a^{a+h} \\&= \dfrac{f(a)-f(a-h)}{2h}\left((a+h-a)^2-(a-a)^2 \right)+f(a)(a+h-a) \\&= \dfrac{f(a)-f(a-h)}{2h} h^2 + hf(a) \\&= h \dfrac{3f(a)-f(a-h)}{2}. \end{aligned}\)

On pose \(a=t_n\) et \(a+h=t_{n+1}\) d’où la formule de quadrature \( \int_{t_{n}}^{t_{n+1}}\!\!\!\!f(t)\mathrm{d}t\approx(t_{n+1}-t_{n}) \frac{3f(t_n)-f(2t_{n}-t_{n+1})}{2}=h \frac{3f(t_n)-f(t_{n-1})}{2}. \) En utilisant la formule de quadrature pour l’intégration de l’EDO \(y'(t)=\varphi(t,y(t))\) entre \(t_n\) et \(t_{n+1}\) on obtient \( y(t_{n+1})=y(t_n)+\int_{t_n}^{t_{n+1}} \varphi(t,y(t))\mathrm{d}t\approx h \dfrac{3\varphi(t_{n},y(t_{n}))-\varphi(t_{n-1},y(t_{n-1}))}{2}. \) Si on note \(u_n\) une approximation de la solution \(y\) au temps \(t_n\), on obtient le schéma à deux pas implicite suivant: \( \begin{cases} u_0=y(t_0)=y_0,\\ u_1\text{ à définir},\\ u_{n+1}=u_n+ h \dfrac{3\varphi(t_{n},u_{n})-\varphi(t_{n-1},u_{n-1})}{2} & n=1,2,\dots N-1 \end{cases} \) On peut utiliser une prédiction d’Euler explicite pour initialiser \(u_1\): \( \begin{cases} u_0=y(t_0)=y_0,\\ u_1=u_0+h\varphi(t_0,u_0),\\ u_{n+1}=u_n+ h \dfrac{3\varphi(t_{n},u_{n})-\varphi(t_{n-1},u_{n-1})}{2} & n=1,2,\dots N-1 \end{cases} \)

Exercice : Construction d’un schéma

Considérons le problème de Cauchy

\(\begin{cases} y'(t) = \varphi(t,y(t)), &\forall t \in I=[t_0,T],\\ y(t_0) = y_0, \end{cases}\)

avec \(y_0\) une valeur donnée et supposons que l’on ait montré l’existence et l’unicité d’une solution \(y\) pour \(t\in I\).

On subdivise l’intervalle \(I=[t_0;T]\), avec \(T<+\infty\), en \(N\) intervalles \([t_{n};t_{n+1}]\) de largeur \(h=(T-t_0)/N\) avec \(t_n=t_0+nh\) pour \(n=0,\dots,N\).
Pour chaque nœud \(t_n\), on note \(y_n=y(t_n)\) et on cherche la valeur inconnue \(u_n\) qui approche la valeur exacte \(y_n\)

Si nous intégrons l’EDO \(y'(t)=\varphi(t,y(t))\) entre \(t_{n-2}\) et \(t_{n+1}\) nous obtenons \( y_{n+1}-y_{n-2}=\int_{t_{n-2}}^{t_{n+1}} \varphi(t,y(t))dt. \) Écrire le schéma obtenu en approchant l’intégrale \(\int_{t_{n-2}}^{t_{n+1}} \varphi(t,y(t))\mathrm{d}t\) par l’intégrale \(\int_{t_{n-2}}^{t_{n+1}} p(t) \mathrm{d}t\) où \(p\) est le polynôme interpolant \(\varphi\) en \(t_{n}\) et \(t_{n-1}\) (attention à bien initialiser la suite définie par récurrence pour qu’on puisse effectivement calculer tous les termes). Le schéma obtenu est-il explicite ?

Correction

Le polynôme interpolant \(\varphi\) en \(t_{n}\) et \(t_{n-1}\) a équation

\( p(t) = \frac{\varphi(t_n,y_n)-\varphi(t_{n-1},y_{n-1})}{t_n-t_{n-1}}(t-t_n)+ \varphi(t_n,y_n). \)

On intègre ce polynôme entre \(t_{n-2}\) et \(t_{n+1}\):

\(\begin{align*} \int_{t_{n-2}}^{t_{n+1}} p(t) \mathrm{d}t &= \frac{\varphi(t_n,y_n)-\varphi(t_{n-1},y_{n-1})}{h}\left[\frac{(t-t_n)^2}{2}\right]_{t_{n-2}}^{t_{n+1}}+ \varphi(t_n,y_n)\left[t\right]_{t_{n-2}}^{t_{n+1}} \\ &=\frac{\varphi(t_n,y_n)-\varphi(t_{n-1},y_{n-1})}{2h}\left[ (t_{n+1}-t_n)^2-(t_{n-2}-t_n)^2 \right]+ \varphi(t_n,y_n)\left[t_{n+1}-t_{n-2}\right] \\ &=\frac{3}{2}h\left(\varphi(t_n,y_n)+\varphi(t_{n-1},y_{n-1})\right). \end{align*}\)

On obtient le schéma explicite

\(\begin{cases} u_0=y_0,\\ u_1=u_0+h\varphi(t_0,u_0), \\ u_{n+1}=u_{n-2}+\frac{3}{2}h\left(\varphi(t_n,u_n)+\varphi(t_{n-1},u_{n-1})\right). \end{cases}\)

import sympy as symb
symb.init_printing()

h = symb.Symbol('h')
t_n = symb.Symbol('t_{n}')
t_nm1 = symb.Symbol('t_{n-1}')
t_np1 = symb.Symbol('t_{n+1}')
t_nm2 = symb.Symbol('t_{n-2}')
f_n = symb.Symbol('f(t_{n})')
f_nm1 = symb.Symbol('f(t_{n-1})')

p = symb.interpolate([(t_n,f_n),(t_nm1,f_nm1)], t).factor(t)
display(symb.Eq(symb.Symbol('p'), p))
q = (symb.integrate(p,(t,t_nm2,t_np1)).subs(t_np1,t_n+h).subs(t_nm1,t_n-h).subs(t_nm2,t_n-2*h)).simplify()
display(symb.Eq( symb.Integral(symb.Symbol('p(t)'), (t, t_nm2,t_np1)) , q ))

\(\displaystyle p = \frac{- f(t_{n-1}) t_{n} + f(t_{n}) t_{n-1} + t \left(f(t_{n-1}) - f(t_{n})\right)}{t_{n-1} - t_{n}}\)

\(\displaystyle \int\limits_{t_{n-2}}^{t_{n+1}} p(t)\, dt = \frac{3 h \left(f(t_{n-1}) + f(t_{n})\right)}{2}\)

Exercice : construction d’un schéma

Considérons le problème de Cauchy:

\( \begin{cases} y'(t) = \varphi(t,y(t)), &\forall t \in [t_0,T],\\ y(t_0) = y_0. \end{cases} \)

Supposons que l’on ait montré l’existence d’une unique solution \(y\).

On subdivise l’intervalle \([t_0,T]\) en \(N\) intervalles de longueur \(h=(T-t_0)/N=t_{n+1}-t_n\). Pour chaque nœud \(t_n=t_0 + nh\) (\(1 \le n \le N\)) on cherche la valeur inconnue \(u_n\) qui approche \(y(t_n)\). L’ensemble de \(N+1\) valeurs \(\{u_0 = y_0, u_1,\dots, u_{N} \}\) représente la solution numérique.

Dans cette exercice on va construire des nouveaux schémas numériques basés sur l’intégration de l’EDO \(y'(t)=\varphi(t,y(t))\) entre \(t_n\) et \(t_{n+2}\): \( y(t_{n+2})=y(t_n)+\int_{t_n}^{t_{n+2}} \varphi(t,y(t))\mathrm{d}t. \)

1. En utilisant la formule de quadrature du point milieu pour approcher le membre de droite écrire un schéma numérique explicite permettant de calculer \(u_{n+2}\) à partir de \(u_{n+1}\) et \(u_{n}\). Notons que ce schéma a besoin de deux valeurs initiales ; on posera alors \(u_0=y_0\) et \(u_{1}\) sera approché par une prédiction d’Euler progressive.
2. En utilisant la formule de quadrature de Cavalieri-Simpson pour approcher le membre de droite écrire un schéma numérique implicite permettant de calculer \(u_{n+2}\) à partir de \(u_{n+1}\) et \(u_{n}\). Notons que ce schéma a besoin de deux valeurs initiales ; on posera alors \(u_0=y_0\) et \(u_{1}\) sera approché par une prédiction d’Euler progressive.
3. Proposer une modification du schéma au point précédent pour qu’il devienne explicite.

Rappels: une formule de quadrature interpolatoire approche l’intégrale d’une fonction \(f\) par l’intègrale dun polynôme qui interpole \(f\) en des points donnés: \(\int_a^b f(t)\mathrm{d}t\approx\int_a^b p(t)\mathrm{d}t\) + la formule de quadrature du point milieu considère comme points d’interpolation le point \(\frac{a+b}{2}\) ainsi \( \int_a^b f(t)\mathrm{d}t\approx\int_a^b f\left(\frac{a+b}{2}\right)\mathrm{d}t = (b-a)f\left(\frac{a+b}{2}\right) \) + la formule de Simpson considère comme points d’interpolation les points \(\left\{a,\frac{a+b}{2},b\right\}\) ainsi \( \int_a^b f(t)\mathrm{d}t\approx \frac{b-a}{6}\left[f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right] \)

Correction premier point

Si on utilise la formule de quadrature du point milieu sur l’intervalle \([t_n;t_{n+2}]\), ie \( \int_{t_n}^{t_{n+2}} \varphi(t,y(t))\mathrm{d}t\approx 2h \varphi\left(t_{n+1},y(t_{n+1})\right) \) on obtient \(\begin{cases} u_0=y(t_0)=y_0,\\ u_1\approx y(t_1)\\ u_{n+2}=u_{n}+2h \varphi(t_{n+1},u_{n+1})& n=0,1,\dots N-2 \end{cases}\) où \(u_{1}\) est une approximation de \(y(t_1)\). Nous pouvons utiliser une prédiction d’Euler progressive pour approcher \(u_1\). Nous avons construit ainsi un nouveau schéma explicite à 2 pas: \(\begin{cases} u_0=y(t_0)=y_0,\\ u_1=u_0+h\varphi(t_{0},u_{0}),\\ u_{n+2}=u_{n}+2h \varphi(t_{n+1},u_{n+1})& n=0,1,\dots N-2 \end{cases}\)

Correction second point

Si on utilise la formule de quadrature de Cavalieri-Simpson sur l’intervalle \([t_n;t_{n+2}]\), ie \( \int_{t_n}^{t_{n+2}} \varphi(t,y(t))\mathrm{d}t\approx \frac{h}{3} \left( \varphi(t_{n},y(t_{n}))+4\varphi(t_{n+1},y(t_{n+1}))+\varphi(t_{n+2},y(t_{n+2})) \right), \) et avec une prédiction d’Euler explicite pour approcher \(u_1\), nous obtenons un nouveau schéma implicite à 3 pas: \(\begin{cases} u_0=y(t_0)=y_0,\\ u_1=u_0+h\varphi(t_{0},u_{0}),\\ u_{n+2}=u_{n}+\frac{h}{3} \left( \varphi(t_{n},u_{n})+4\varphi(t_{n+1},u_{n+1})+\varphi(t_{n+2},u_{n+2}) \right)& n=0,1,\dots N-2 \end{cases}\)

Correction troisième point

Pour éviter le calcul implicite de \(u_{n+2}\), nous pouvons utiliser une technique de type predictor-corrector et remplacer le \(u_{n+2}\) dans le terme \(\varphi(t_{n+2},u_{n+2})\) par exemple par \(\tilde u_{n+2}=u_{n+1}+h \varphi(t_{n+1},u_{n+1})\). Nous avons construit ainsi un nouveau schéma explicite à 3 pas: \(\begin{cases} u_0=y(t_0)=y_0,\\ u_1=u_0+h\varphi(t_{0},u_{0}),\\ u_{n+2}=u_{n}+\frac{h}{3} \left( \varphi(t_{n},u_{n})+4\varphi(t_{n+1},u_{n+1})+\varphi(t_{n+2},u_{n+1}+h \varphi(t_{n+1},u_{n+1}) \right)& n=0,1,\dots N-2 \end{cases}\)