from IPython.core.display import HTML
css_file = './custom.css'
HTML(open(css_file, "r").read())

import sys #only needed to determine Python version number
print('Python version ' + sys.version)

Python version 3.8.5 (default, Jul 28 2020, 12:59:40) 
[GCC 9.3.0]

M62_TD1 : Exercices 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$.

%reset -f
%matplotlib inline

from matplotlib.pylab import *

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

phi = lambda t,y : a*y # EDO
sol_exacte = lambda t : y0*math.exp(a*t) # SOLUTION EXACTE

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

# CALCUL SOLUTION APPROCHEE
tt1 = linspace(t0,tfinal,Nh+1)
uu1 = EE(phi,tt1,y0)

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

# PLOT
plot(tt1,yy,'m-',tt1,uu1,'go-') 	
legend(['Exacte','Euler explicite'],loc='upper left')

# PRINT
h1=tt1[1]-tt1[0]
for n in range(Nh):
	print ('y_%1d=%1.14f u_%1d=%1.14f (1+h)^%1d=%1.14f' %(n,yy[n],n,uu1[n],n,(1.+h1)**n))

y_0=1.00000000000000 u_0=1.00000000000000 (1+h)^0=1.00000000000000
y_1=1.49182469764127 u_1=1.40000000000000 (1+h)^1=1.40000000000000
y_2=2.22554092849247 u_2=1.96000000000000 (1+h)^2=1.96000000000000
y_3=3.32011692273655 u_3=2.74400000000000 (1+h)^3=2.74400000000000
y_4=4.95303242439511 u_4=3.84160000000000 (1+h)^4=3.84160000000000

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

# MAIN
y_2=sol_exacte(tfinal)
print ('\nSolution exacte en t=2:')
print ('y(2) =',y_2)
print ('\nSolution approchee en t=2 avec N points:')
print ('N\t sol_approx\t\t\t |erreur|')
for i in range(3,15):
    N=2**i
    tt = linspace(t0,tfinal,N+1)
    uu = EE(phi,tt,y0)
    print (N,'\t',uu[-1], '\t',abs(uu[-1]-y_2))

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

Solution approchee en t=2 avec N points:
N	 sol_approx			 |erreur|
8 	 5.9604644775390625 	 1.428591621391588
16 	 6.583250172027423 	 0.8058059269032274
32 	 6.958666757218805 	 0.4303893417118454
64 	 7.166276152788222 	 0.22277994614242846
128 	 7.275669793128417 	 0.11338630580223352
256 	 7.3318505987410365 	 0.0572055001896139
512 	 7.3603235532692795 	 0.028732545661370956
1024 	 7.374657160341845 	 0.014398938588805699
2048 	 7.381848435880501 	 0.007207663050149193
4096 	 7.385450215539008 	 0.0036058833916428057
8192 	 7.38725264383889 	 0.0018034550917604975
16384 	 7.388154242982074 	 0.0009018559485767241

Vérifions la convergence linéaire:

error=[]
H=[]
N=arange(200,1000,100)
for n in N:
    tt=linspace(t0,tfinal,n+1)
    H.append(tt[1]-tt[0])
    yy=[sol_exacte(t) for t in tt]
    uu = EE(phi,tt,y0)
    error_vec=[abs(uu[i]-yy[i]) for i in range(len(uu))]
    error.append(max(error_vec))

# ESTIMATION ORDRES DE CONVERGENCE par regression
print ('Pente par regression lineaire', polyfit(log(H),log(error), 1)[0])

# PLOT ERREUR     
loglog(H,error,'*-');

Pente par regression lineaire 0.9942137370997959

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 obtenue en appliquant la méthode d'\textsc{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]$ pour les valeurs de $y_0$ suivantes: \begin{align*} 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{align*}
   Sur un autre graphique à coté, tracer les solutions correspondantes obtenues par la méthode d'Euler explicite avec $N_h=50$.
   Sur un autre graphique à coté, tracer les solutions correspondantes obtenues par la méthode d'Euler explicite avec $N_h=100$.

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 environnementale (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}$$

from matplotlib.pylab import *

def EE(phi,tt,y0):
    h=tt[1]-tt[0]
    uu = [y0]
    for i in range(len(tt)-1):
        uu.append(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
        return 0
    elif abs(y0-a/b)<=1.e-10:
        return a/b
    else:
        return 1./((1./y0-b/a)*exp(-a*t)+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]
leg=[]
for y0 in yy_0:
	leg.append('y_0='+str(y0))

# MAIN
figure(1, figsize=(15, 5))
subplot(1,3,1)
tt=linspace(t0,tfinal,25)
for i in range(len(yy_0)):
	uu=EE(phi,tt,yy_0[i])
	plot(tt,uu,'.-')
axis([t0,tfinal,0,2.])
title('25 noeuds')
grid()
legend(leg)

subplot(1,3,2)
tt=linspace(t0,tfinal,100)
for i in range(len(yy_0)):
	uu=EE(phi,tt,yy_0[i])
	plot(tt,uu,'.-')
title('100 noeuds')
axis([t0,tfinal,0,2.])
grid()
legend(leg)

subplot(1,3,3)
for i in range(len(yy_0)):
	yy=[sol_exacte(t,yy_0[i]) for t in tt]
	plot(tt,yy,'-')
axis([t0,tfinal,0,2.])
title('Exacte')
grid()
legend(leg);

 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}$.

%matplotlib inline
from matplotlib.pylab import *
y0=1
exacte = lambda t : y0*exp(-10*t)
tt=linspace(0,3,101)
yy=[exacte(t) for t in tt]
plot(tt,yy);

Correction

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}. $$

Correction

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

Correction

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

%matplotlib inline
from matplotlib.pylab import *

N=6
tt=linspace(0,1,N+1)

y0=1
exacte = lambda t : y0*exp(-10*t)
yy=[exacte(t) for t in tt]

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

plot(tt,yy,label=("Exacte"))
plot(tt,schema(phi,tt,y0),"-.",label=("Schéma"))
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}$$

1. 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}, $$ pour $h>0$ fixé. 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}. $$

Correction

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

%matplotlib inline
from matplotlib.pylab import *

y0=1

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

N=10
tt=linspace(0,10,N+1)

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*sqrt(uu[i])))
    return uu

yy=[exacte(t) for t in tt]
plot(tt,yy,'r-',label=('Exacte'))
plot(tt,schema(phi,tt,y0),'b*',label=("Schéma"))
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} $$

y0=1
exacte = lambda t : y0*(4*t*y0**4+1)**(-0.25)
tt=linspace(0,10,101)
yy=[exacte(t) for t in tt]
plot(tt,yy);

Correction

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}. $$

Correction

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.

N=50
tt=linspace(0,10,N+1)

y0=1
exacte = lambda t : y0*(4*t*y0**4+1)**(-0.25)
yy=[exacte(t) for t in tt]

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

plot(tt,yy,label=("Exacte"),lw=2)
plot(tt,schema(phi,tt,y0),'o',label=("Schéma"))
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.

Correction

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}. $$

Correction

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

Correction

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

N=100
tt=linspace(0,10,N+1)

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

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

plot(tt,schema(phi,tt,y0),'o',label=("Schéma"))
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}$$ Montrer que si $0<I_0<A$ alors $0<I(t)<A$ pour tout $t>0$.
   Montrer que si $0<I_0<A$ alors $I(t)$ est croissante sur $\mathbb{R}^+$.
2. 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é.
3. 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}^+$.

Correction

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.

Correction

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{align*} % 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{align*}

La valeur numérique de la constante d'intégration $D$ est obtenue grâce à la CI: \begin{align*} % I(0)&=\frac{A}{De^{0}+1}\\ D&=\frac{A-I_0}{I_0} \end{align*}

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

A=5000
I0=160
k=log(363/38)/35000

D=(A-I0)/I0

exacte = lambda t : A/(D*exp(-A*k*t)+1)
h=1
tt=arange(0,30,h)
yy=[exacte(t) for t in tt]

phi = lambda t,I : k*I*(A-I) 

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

plot(tt,yy,label=("Exacte"))
plot(tt,schema(phi,tt,y0),label=("Schéma"))
legend();

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}$$

1. Écrire le problème comme un système de 2 EDO d'ordre 1 et en calculer la solution.
2. Écrire le schéma le schéma d'Euler explicite pour ce problème.
3. Pour quelles valeurs de $h$ le schéma est A-stable?

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!

%reset -f
import sympy as symb
symb.init_printing()

symb.var('t')
M = symb.Matrix(((0,-1), (1000,1001)))
M

$\displaystyle \left[\begin{matrix}0 & -1\\1000 & 1001\end{matrix}\right]$

p=M.charpoly(t)
symb.factor(p)

$\displaystyle \left(t - 1000\right) \left(t - 1\right)$

P, D = M.diagonalize()
P,D

$\displaystyle \left( \left[\begin{matrix}-1 & -1\\1 & 1000\end{matrix}\right], \ \left[\begin{matrix}1 & 0\\0 & 1000\end{matrix}\right]\right)$

symb.simplify(P*D*P**-1)

$\displaystyle \left[\begin{matrix}0 & -1\\1000 & 1001\end{matrix}\right]$

Correction

Soit $h>0$ fixé et $t_n=nh$ pour tout $n\in\mathbb{Z}$. Le schéma d'Euler explicite pour le système donné construit la suite $$ \begin{aligned} \mathbf{u}_{n+1} &= \mathbf{u}_n -h \begin{pmatrix} 0&-1\\ 1000&1001 \end{pmatrix} \mathbf{u}_n \\ &= \begin{pmatrix} 1 & h\\ -1000h & 1-1001h \end{pmatrix} \mathbf{u}_n \\ &= \begin{pmatrix} 1 & h\\ -1000h & 1-1001h \end{pmatrix}^{n+1} \mathbf{u}_0,\\ &= \begin{pmatrix} -1 & -1\\ 1000 & 1 \end{pmatrix} \begin{pmatrix} 1-1000h & 0\\ 0 & 1-h \end{pmatrix}^{n+1} \begin{pmatrix} -1 & -1\\ 1000 & 1 \end{pmatrix}^{-1} \mathbf{u}_0\\ &= \begin{pmatrix} -1 & -1\\ 1000 & 1 \end{pmatrix} \begin{pmatrix} (1-1000h)^{n+1} & 0\\ 0 & (1-h)^{n+1} \end{pmatrix} \begin{pmatrix} -1 & -1\\ 1000 & 1 \end{pmatrix}^{-1} \mathbf{u}_0, \qquad \forall n \in \mathbb{N}. \end{aligned} $$

symb.var('h')
M = symb.Matrix(((1,h), (-1000*h,1-1001*h)))
M

$\displaystyle \left[\begin{matrix}1 & h\\- 1000 h & 1 - 1001 h\end{matrix}\right]$

P, D = M.diagonalize()
P , D

$\displaystyle \left( \left[\begin{matrix}-1 & -1\\1000 & 1\end{matrix}\right], \ \left[\begin{matrix}1 - 1000 h & 0\\0 & 1 - h\end{matrix}\right]\right)$

Correction

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}$.

Exercice : A-stabilité du $\vartheta$-schéma

Pour calculer la solution approchée 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 suivante: $$\begin{cases} u_0=y_0,\\ u_{n+1}=u_{n}+h\left(\vartheta\varphi(t_{n+1},u_{n+1})+(1-\vartheta)\varphi(t_{n},u_{n})\right). \end{cases}$$ avec $\vartheta \in [0;1]$.

Soit $\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.

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-2. 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} $$

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:

%reset -f
%matplotlib inline
%autosave 300

from matplotlib.pylab import *

H=[0.1,0.05,0.025]
err=[]
y=0.460857
err.append(abs(y-0.459861))
err.append(abs(y-0.460633))
err.append(abs(y-0.460804))

print ('Ordre\t %1.2f' %(polyfit(log(H),log(err), 1)[0]))

subplot(1,2,1)
loglog(H,err, 'r-o');
grid()
subplot(1,2,2)
plot(log(H),log(err), 'r-o');

Autosaving every 300 seconds
Ordre	 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 \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 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()

symb.var('f_0,f_1,f_m1,t')

p=symb.interpolate([(1,f_1),(0,f_0),(-1,f_m1)], t).factor(t)
p

$\displaystyle - \frac{- 2 f_{0} + t^{2} \left(2 f_{0} - f_{1} - f_{m1}\right) + t \left(- f_{1} + f_{m1}\right)}{2}$

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}. $$

q=(symb.integrate(p,(t,0,1)).simplify())
q

$\displaystyle \frac{2 f_{0}}{3} + \frac{5 f_{1}}{12} - \frac{f_{m1}}{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()

symb.var('a,h,f_a,f_amh,t')

p=symb.interpolate([(a,f_a),(a-h,f_amh)], t).factor(t)
p

$\displaystyle \frac{- a f_{a} + a f_{amh} + f_{a} h + t \left(f_{a} - f_{amh}\right)}{h}$

On en déduit la méthode de quadrature \begin{align*} \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{align*}

q=(symb.integrate(p,(t,a,a+h)).simplify())
q

$\displaystyle \frac{h \left(3 f_{a} - f_{amh}\right)}{2}$

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

trouver une fonction $y \colon I\subset \mathbb{R} \to \mathbb{R}$ définie sur un intervalle $I=[t_0,T]$ telle que $$\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}$$