TD1 - Exercices sur les schémas à un pas – M64 – Data Science & Scientific Computing 6 : Schémas numériques pour les problèmes de Cauchy

TD1 - Exercices sur les schémas à un pas

Auteur·rice

Gloria Faccanoni

Date de publication

20 mars 2026

Code 
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import fsolve
from scipy.integrate import solve_ivp, odeint
import sympy as sp
from IPython.display import display, Math, Markdown
# plt.rcParams.update({'font.size': 12})
plt.rcdefaults() # Réinitialisation des paramètres par défaut de Matplotlib
# Mise à jour des paramètres de police
# plt.rcParams.update({'font.size': 12})
# plt.rcParams.update({'figure.max_open_warning': 0})

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\).

Code 
def EE(phi, tt, y0):
    h = tt[1] - tt[0]
    uu = np.empty_like(tt)
    uu[0] = y0
    for n in range(len(tt) - 1):
        uu[n + 1] = uu[n] + h * phi(tt[n], uu[n])
    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)  

# PLOT
plt.plot(tt, yy, 'm-', label='Exacte')
plt.plot(tt, uu, 'go-', label='Euler explicite')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
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). \)

Code 
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:

Code 
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=(8, 5))
plt.title(f'Ordre = {pente}')
plt.loglog(H, error, '*-');

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}\)

Code 
def EE(phi, tt, y0):
    h = tt[1] - tt[0]
    uu = np.empty_like(tt)
    uu[0] = y0
    for n in range(len(tt) - 1):
        uu[n + 1] = uu[n] + h * phi(tt[n], uu[n])
    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'])

# Version directement vectorisée sans np.vectorize
sol_exacte_vectorized = lambda y, y0 : np.where(  abs(y0-0)<=1.e-10, 0 , np.where( abs(y0 - a / b) <= 1.e-10, a / b, 1 / ((1 / y0 - b / a) * np.exp(-a * y) + b / a) )  )

# 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()

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\).

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}\).

Code 
# 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}. \)

Code 
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 n in range(len(tt)-1):
        uu[n+1] = (1 - 5*h) / (1 + 5*h) * uu[n]
    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();

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\).

Code 
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 n in range(len(tt)-1):
        uu.append(2 * (uu[n]) ** (3/2) / (h + 2 * np.sqrt(uu[n])))
    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();

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} \)

Code 
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.

Code 
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 n in range(len(tt) - 1):
        uu[n + 1] = uu[n] / (1 + h * uu[n] ** 4)
    return uu

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

plt.figure(figsize=(8, 5))
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();

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\).

Code 
t0 = 0 
y0 = 1

# SOLUTION EXACTE POUR y(0) = 1

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

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

solgen = solgenLIST[1]
display(solgen)

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

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

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

func = sp.lambdify(t,solpar.rhs,'numpy')
[Eq(y(t), -acos((-C1 - exp(2*t))/(C1 - exp(2*t))) + 2*pi),
 Eq(y(t), acos((-C1 - exp(2*t))/(C1 - exp(2*t))))]

\(\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)}\)

{C1: (-1 + cos(1))/(cos(1) + 1)}

\(\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)}\)

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

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
Code 
# COMPARAISONS
N = 100
tt = np.linspace(t0,10,N+1)

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

plt.figure(figsize=(8,5))
plt.plot(tt, yy, label="Exacte")
plt.plot(tt, uu, 'o', label="Schéma")
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left');

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\):

Code 
# 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 n in range(len(tt) - 1):
        uu[n + 1] = (1 + k * A * h) / (1 + k * uu[n] * h) * uu[n]
    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();

\( 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}\).

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:

Code 
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, 5))
# 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