2022 CR

• ⏱ 2H30
• ✍ Deposer votre notebook .ipynb sur Moodle (dépôt de devoir)

Exercice 1 : Méthode d'Euler exponentielle et problèmes raides

Soient $b,g>0$ deux constantes et considérons le problème de Cauchy $$\begin{cases}y'(t)=-by(t)+g,& t\in[0;1]\\y(0)=1+10^{-8}.\end{cases}$$

Q1 [2 points] Solution exacte
Calculer et tracer le graphe de la solution exacte $t\mapsto y(t)$ d'abord pour $b=g=10$ puis pour $b=g=100$.

Il s'agit d'une edo linéaire d'ordre 1 dont la solution est $$ y(t)=\left( y_0-\frac{g}{b} \right)e^{-bt}+\frac{g}{b} \xrightarrow[t\to+\infty]{}y(t)=\frac{g}{b} \text{ pour tout } y(0)=y_0\in\mathbb{R} $$ Le problème est numériquement bien posé, c'est-à-dire que l'erreur sur la solution sera au pire de l'ordre de l'erreur sur la condition initiale. Cependant, lorsque $b$ devient de plus en plus grand, la solution devient de plus en plus raide: il s'agit d'un problème stiff.

%reset -f
%matplotlib inline
from matplotlib.pylab import *

# ======================================================
# solution exacte
# ======================================================

exacte = lambda t,b,g,y0 : (y0-g/b)*exp(-b*t)+g/b 

# ======================================================
# variables globales
# ======================================================
t0 = 0
tfinal = 1
y0 = 1+1.e-8

# ======================================================
# affichage
# ======================================================
tt = linspace(t0,tfinal,101)
figure(figsize=(10,7))
for b in [10,100] :
    plot( tt , [exacte(t,b,b,y0) for t in tt] , label=rf'$b=g={b}$');
legend()
grid()

Q1 [4 points] Méthode d'Euler exponentielle La méthode d'Euler exponentielle appliquée à ce problème de Cauchy s'écrit comme suit: $$\text{(EEx) }\begin{cases}u_0=y_0,\\u_{n+1}=e^{-bh}u_n+ghf(-bh), & n=0,\dots, N-1\end{cases}$$ où $f(z)=\frac{e^z-1}{z}$.

1. Vérifier que cette méthode génère une suite arithmético-géométrique et montrer analytiquement qu'elle calcule la solution exacte discrète.
2. Soient $b=g=100$ et $y_0=1+10^{-8}$. Tracer la solution exacte et les solutions approchées obtenues avec $N=10$, $N=50$, $N=100$.

La méthode d'Euler exponentielle donne la suite définie par récurrence $$ \begin{cases} u_0=y_0,\\ u_{n+1}=Au_n+B \end{cases} \qquad\text{avec } A=\exp(-bh), \ B=(\exp(-bh)-1)\frac{g}{b}=(A-1)\frac{g}{b} $$ Donc $$ u_{n+1}=Au_n+(A-1)\frac{g}{b}=A\left(u_n+\frac{g}{b}\right)-\frac{g}{b} $$ soit encore $$ u_{n+1}+\frac{g}{b} = A\left(u_n+\frac{g}{b}\right). $$ Si on note $w_n=u_n+\frac{g}{b}$, on a $$ w_{n+1}=Aw_n=A^{n+1}w_0 $$ ainsi $$ u_n-\frac{g}{b}=(\exp(-hb))^n\left( y_0 -\frac{g}{b}\right). $$ Comparons la solution exacte et la solution approchée: \begin{align*} y_n&=\frac{g}{b}+\left( y_0 -\frac{g}{b}\right)e^{-bt_n}\\ u_n&=\frac{g}{b}+\left( y_0 -\frac{g}{b}\right)e^{-bhn} \end{align*} Comme $t_n=t_0+nh=nh$, la solution approchée est exacte.

b = 100
g = b

# ======================================================
# SCHÉMA
# ======================================================

def EEx(tt,y0):
    uu = [y0]
    h = tt[1]-tt[0]
    for i in range(len(tt)-1):
        K1 = (exp(-b*h)-1)/(-b*h)
        uu.append( exp(-b*h)*uu[i]+g*h*K1 )
    return uu


# ======================================================
# CALCUL SOLUTION APPROCHEE
# ======================================================
figure(1,figsize=(30,7))
idx=1
for N in [10,50,100]:
    subplot(1,3,idx)
    tt = linspace(t0,tfinal,N+1)
    yy = [exacte(t,b,b,1+1.e-8)  for t in tt]
    uu = EEx(tt,y0)
    plot(tt,yy,'-',label='Exacte')
    plot(tt,uu,'-D',label=f'{N=}')    
    legend()
    grid()
    idx+=1

Exercice 2 : étude d'un schéma Runge-Kutta

Considérons le schéma dont la matrice de Butcher est donnée ci-dessous: $$ \begin{array}{c|cccc} \dfrac{1}{2} & \dfrac{1}{2} & 0 & 0 & 0\\[1em] \dfrac{2}{3} & \dfrac{1}{6} & \dfrac{1}{2}& 0 & 0\\[1em] \dfrac{1}{2} & -\dfrac{1}{2} & \dfrac{1}{2}& \color{red}\vartheta & 0\\[1em] 1 & \dfrac{3}{2} & -\dfrac{3}{2}& \dfrac{1}{2} & \dfrac{1}{2}\\[1em] \hline & \dfrac{3}{2} & \color{red}\gamma & \dfrac{1}{2} & \dfrac{1}{2} \end{array} $$

Q1 [1 points] Calculer $\gamma$ et $\vartheta$ pour que le schéma soit consistant. Dans la suite on posera $\gamma$ et $\vartheta$ égales à ces valeurs.

C'est un schéma semi-implicite à $4$ étages ($s=4$).

Pour qu'il soit consistant il faut que $\begin{cases} \displaystyle\sum_{j=1}^s b_{i}=1& \\ \displaystyle c_i=\sum_{j=1}^s a_{ij}&i=1,\dots,s \end{cases}$

• Comme $\displaystyle\sum_{j=1}^s b_{i}=\dfrac{3}{2} + \gamma + \dfrac{1}{2} + \dfrac{1}{2}=\dfrac{5}{2} + \gamma$, il faut poser $\gamma=-\dfrac{3}{2}$.
• La deuxième condition est satisfaite ssi $\dfrac{1}{2} = -\dfrac{1}{2} +\dfrac{1}{2}+ \vartheta $, il faut donc poser $\vartheta=\dfrac{1}{2} $.

Q2 [5 points] Étudier empiriquement l'ordre de convergence sur le problème de Cauchy $$\begin{cases} y'(t) = -\dfrac{y(t)}{2}\left(1-\dfrac{y(t)}{3}\right), &\forall t \in I=[0,10],\\ y(0) = 1, \end{cases} $$

%reset -f
%matplotlib inline

from matplotlib.pylab import *
from scipy.optimize import fsolve

import sympy 
sympy.init_printing()

from IPython.display import display, Math
prlat= lambda *args: display(Math(''.join( sympy.latex(arg) for arg in args )))

# variables globales
t0     = 0
tfinal = 10
y0     = 1

phi = lambda t,y: -y/2*(1-y/3) #-y-t**2

##############################################
# solution exacte
##############################################
t   = sympy.Symbol('t')
y   = sympy.Function('y')
edo = sympy.Eq( sympy.diff(y(t),t) , phi(t,y(t)) )
solgenLIST = sympy.dsolve(edo)
display(solgenLIST)
solgen=solgenLIST[0]
display(solgen)
consts = sympy.solve( sympy.Eq( y0, solgen.rhs.subs(t,t0)) , dict=True)[0]
display(consts)
solpar = solgen.subs(consts).simplify()
display(solpar)

sol_exacte = sympy.lambdify(t,solpar.rhs,'numpy')

##############################################
# solution approchée
##############################################
def RK(tt,phi,y0):
    uu = [y0]
    h = tt[1]-tt[0]
    for i in range(len(tt)-1):
        sys = lambda X : [ -X[0]+phi( tt[i]+h/2   ,uu[i]+h/2*X[0]                      ) , 
                           -X[1]+phi( tt[i]+h*2/3 ,uu[i]+h/6*(X[0]+3*X[1])             ), 
                           -X[2]+phi( tt[i]+h/2   ,uu[i]+h/2*(-X[0]+X[1]+X[2])         ), 
                           -X[3]+phi( tt[i+1]     ,uu[i]+h/2*(3*X[0]-3*X[1]+X[2]+X[3]) )
                           ]
        K_start = phi(tt[i],uu[i])
        K1, K2, K3, K4 = fsolve( sys , [K_start,K_start,K_start,K_start] )
        uu.append( uu[i]+h/2*(3*K1-3*K2+K3+K4) )
    return uu


##############################################
# ordre
##############################################
H = []
err = []
N = 10

figure(figsize=(16,5))
ax1 = subplot(1,2,1)
ax2 = subplot(1,2,2)

for k in range(6):
    N += 20
    tt = linspace(t0, tfinal, N + 1)
    H.append( tt[1] - tt[0] )
    yy = sol_exacte(tt)
    uu = RK(tt,phi,y0)
    err.append( norm(uu-yy,inf) )
    ax1.plot(tt,uu,label=f'Approchée avec N={N}')

ax1.plot(tt,yy,label='Exacte')
xlabel('$t$')
ylabel('$y$')
ax1.grid(True)
ax1.legend()

ax2.loglog( H, err, 'r-o')
xlabel('$h$')
ylabel('$e$')
title(f'RK : ordre = {polyfit(log(H),log(err),1)[0]:1.2f}')
ax2.grid(True)

$\displaystyle \left[ y{\left(t \right)} = \frac{3 \left(e^{C_{1}} - \sqrt{e^{C_{1} + t}}\right)}{e^{C_{1}} - e^{t}}, \ y{\left(t \right)} = \frac{3 \left(e^{C_{1}} + \sqrt{e^{C_{1} + t}}\right)}{e^{C_{1}} - e^{t}}\right]$

$\displaystyle y{\left(t \right)} = \frac{3 \left(e^{C_{1}} - \sqrt{e^{C_{1} + t}}\right)}{e^{C_{1}} - e^{t}}$

$\displaystyle \left\{ C_{1} : - \log{\left(4 \right)}\right\}$

$\displaystyle y{\left(t \right)} = \frac{3 \left(2 \sqrt{e^{t}} - 1\right)}{4 e^{t} - 1}$

Q3 [3 points] Étudier théoriquement l'ordre de convergence du schéma.

Soit $\omega$ l'ordre de la méthode.

• C'est un schéma semi-implicite à $s=4$ étages donc $\omega\le2s=8$

• Pour $\vartheta=\frac{1}{2}$ et $\gamma=\frac{3}{2}$, $\omega\ge1$.

• Si, de plus, $\displaystyle\sum_{j=1}^s b_j c_j=\frac{1}{2}$ alors $\omega\ge2$

• Si, de plus,
  $\begin{cases} \displaystyle\sum_{j=1}^s b_j c_j^2=\frac{1}{3}\\ \displaystyle\sum_{i=1}^s\sum_{j=1}^s b_i a_{ij} c_j=\frac{1}{6} \end{cases}$ alors $\omega\ge3$

• Si, de plus, $\begin{cases} \displaystyle\sum_{j=1}^s b_j c_j^3=\frac{1}{4}& \\ \displaystyle\sum_{i=1}^s\sum_{j=1}^s b_i c_i a_{ij} c_j=\frac{1}{8} \\ \displaystyle\sum_{i=1}^s\sum_{j=1}^s b_i a_{ij} c_j^2=\frac{1}{12} \\ \displaystyle\sum_{i=1}^s\sum_{j=1}^s\sum_{k=1}^s b_i a_{ij}a_{jk} c_k=\frac{1}{24} \end{cases}$ alors $\omega\ge4$

Calculons donc toutes ces sommes:

%reset -f

import sympy 
sympy.init_printing()

from IPython.display import display, Math
prlat= lambda *args: display(Math(''.join( sympy.latex(arg) for arg in args )))

c = [sympy.Rational(1,2), sympy.Rational(2,3), sympy.Rational(1,2), 1]
b = [ sympy.Rational(3,2), -sympy.Rational(3,2), sympy.Rational(1,2), sympy.Rational(1,2)]
A = [ [sympy.Rational(1,2), 0, 0, 0], 
      [sympy.Rational(1,6), sympy.Rational(1,2), 0, 0],
      [-sympy.Rational(1,2), sympy.Rational(1,2), sympy.Rational(1,2), 0],
      [sympy.Rational(3,2), -sympy.Rational(3,2), sympy.Rational(1,2), sympy.Rational(1,2)]]
s=len(c)

# Affichage matrice de Butcher
print(f'{"Matrice de Butcher":_^150}')
But = sympy.Matrix(A)
But=But.col_insert(0,sympy.Matrix(c))   
last=[0]
last.extend(b)
But=But.row_insert(s,sympy.Matrix(last).T)
prlat(sympy.Matrix(But))
#################

print(f'{"Consistance i.e. ordre 1 ssi ils sont tous = 0":_^150}')
prlat("\sum_{j=1}^s b_j-1 =",sum(b)-1)
for i in range(s):
    prlat(f'i={i}','\quad \sum_{j=1}^s a_{ij}-c_i=',(sum(A[i])-c[i]))

print(f'{"Ordre 2  ssi (ordre 1 et = 1/2)":_^150}')
prlat('\sum_{j=1}^s b_j c_j=',sum([b[i]*c[i] for i in range(s)]))

print(f'{"Ordre 3 ssi (ordre 2 et premier = 1/3 et, de plus, deuxième 1/6)":_^150}')
prlat('\sum_{j=1}^s b_j c_j^2=',sum([b[i]*c[i]**2 for i in range(s)]).simplify())
prlat('\sum_{i,j=1}^s b_i a_{ij} c_j=',sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]))
print(f'{"Ordre 4 ssi (ordre 3 et, de plus, premier = 1/4, deuxieme = 1/8, troisieme = 1/12 et quatrième = 1/24)":_^150}')
prlat('\sum_{j=1}^s b_j c_j^3=',sum([b[i]*c[i]**3 for i in range(s)]).simplify())
prlat('\sum_{i,j=1}^s b_i c_i a_{ij} c_j=',sum([b[i]*c[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]))
prlat('\sum_{i,j=1}^s b_i a_{ij} c_j^2=',sum([b[i]*A[i][j]*c[j]**2 for i in range(s) for j in range(s)]))
prlat('\sum_{i,j,k=1}^s b_i a_{ij}a_{jk} c_k=',sum([b[i]*A[i][j]*A[j][k]*c[k] for i in range(s) for j in range(s) for k in range(s)]))

_________________________________________________________________Matrice de Butcher :_________________________________________________________________

$\displaystyle \left[\begin{matrix}\frac{1}{2} & \frac{1}{2} & 0 & 0 & 0\\\frac{2}{3} & \frac{1}{6} & \frac{1}{2} & 0 & 0\\\frac{1}{2} & - \frac{1}{2} & \frac{1}{2} & \frac{1}{2} & 0\\1 & \frac{3}{2} & - \frac{3}{2} & \frac{1}{2} & \frac{1}{2}\\0 & \frac{3}{2} & - \frac{3}{2} & \frac{1}{2} & \frac{1}{2}\end{matrix}\right]$

___________________________________________________Consistance i.e. ordre 1 ssi ils sont tous = 0 :___________________________________________________

$\displaystyle \sum_{j=1}^s b_j-1 =0$

$\displaystyle i=0\quad \sum_{j=1}^s a_{ij}-c_i=0$

$\displaystyle i=1\quad \sum_{j=1}^s a_{ij}-c_i=0$

$\displaystyle i=2\quad \sum_{j=1}^s a_{ij}-c_i=0$

$\displaystyle i=3\quad \sum_{j=1}^s a_{ij}-c_i=0$

___________________________________________________________Ordre 2  ssi (ordre 1 et = 1/2)____________________________________________________________

$\displaystyle \sum_{j=1}^s b_j c_j=\frac{1}{2}$

___________________________________________Ordre 3 ssi (ordre 2 et premier = 1/3 et, de plus, deuxième 1/6)___________________________________________

$\displaystyle \sum_{j=1}^s b_j c_j^2=\frac{1}{3}$

$\displaystyle \sum_{i,j=1}^s b_i a_{ij} c_j=\frac{1}{6}$

________________________Ordre 4 ssi (ordre 3 et, de plus, premier = 1/4, deuxieme = 1/8, troisieme = 1/12 et quatrième = 1/24)________________________

$\displaystyle \sum_{j=1}^s b_j c_j^3=\frac{11}{36}$

$\displaystyle \sum_{i,j=1}^s b_i c_i a_{ij} c_j=\frac{5}{48}$

$\displaystyle \sum_{i,j=1}^s b_i a_{ij} c_j^2=\frac{5}{72}$

$\displaystyle \sum_{i,j,k=1}^s b_i a_{ij}a_{jk} c_k=\frac{1}{48}$

D'après ces résultats, le schéma est d'ordre $3$.

Exercice 3 : étude d'un schéma multipas

On notera $\varphi_k\stackrel{\text{déf}}{=}\varphi(t_k,u_k)$. Considérons les deux méthodes multipas $$ u_{n+1}=\frac{18}{11}u_n-\frac{9}{11}u_{n-1}+\frac{2}{11}u_{n-2}+h\frac{6}{11}\varphi_{n+1} $$ $$ u_{n+1}=\frac{18}{11}u_n-\frac{9}{11}u_{n-1}+\frac{3}{11}u_{n-2}+h\frac{6}{11}\varphi_{n+1} $$

Q1 [1 points] Choisir la méthode consistente (justifier la réponse). Dans la suite on n'étudiera que cette méthode.

• On a $p=2$: c'est des méthodes à $q=p+1=3$ pas.
• Les méthodes sont implicites car $b_{-1}=\frac{6}{11}\neq0$.

Pour que la méthode soit consistante il faut que $$\begin{cases} \displaystyle\sum_{j=0}^p a_j=a_0+a_1=1, \\ \displaystyle\sum_{j=0}^p (-j)a_j+\sum_{j=-1}^p b_j=1. \end{cases}$$ Seul la première des méthodes est consistante car $$ \frac{18}{11}-\frac{9}{11}+\frac{2}{11}=1 $$ tandis que $$ \frac{18}{11}-\frac{9}{11}+\frac{3}{11}=\frac{12}{11}\neq1. $$

Q2 [2 points] Étudier la zéro-stabilité?

Le premier polynôme caractéristique est $$ \varrho(r)= r^{p+1}-\sum_{j=0}^p a_jr^{p-j}= r^3-\frac{18}{11}r^2+\frac{9}{11}r-\frac{2}{11}= \frac{1}{11}(r-1)(11r^2-7r+2). $$ Calculons ses racines avec sympy.

%reset -f
%matplotlib inline

from matplotlib.pylab import *
from scipy.optimize import fsolve

import sympy 
sympy.init_printing()

from IPython.display import display, Math
prlat= lambda *args: display(Math(''.join( sympy.latex(arg) for arg in args )))

sympy.var('r')
rho = 11*r**3-18*r**2+9*r-2
sol=sympy.solve(rho)
prlat(sol)
modules=[abs(s) for s in sol]
prlat(modules)

$\displaystyle \left[ 1, \ \frac{7}{22} - \frac{\sqrt{39} i}{22}, \ \frac{7}{22} + \frac{\sqrt{39} i}{22}\right]$

$\displaystyle \left[ 1, \ \frac{\sqrt{22}}{11}, \ \frac{\sqrt{22}}{11}\right]$

Les racines sont $$ r_0=1,\quad r_{1,2}=\frac{7\pm\sqrt{-39}}{22}. $$ La méthode est donc zéro-stable car $$ \begin{cases} |r_j|\le1 \quad\text{pour tout }j=0,1,2 \\ \varrho'(r_j)\neq0 \text{ si } |r_j|=1. \end{cases} $$

Q3 [3 points] Étudier théoriquement l'ordre.

La première barrière de Dahlquist affirme qu'un schéma implicite à $q=3$ pas consistante et zéro-stable peut être au mieux d'ordre $\omega=q+1=4$.

Pour que la méthode soit

• au moins d'ordre 2 il faut qu'elle soit consistante et que $\displaystyle\sum_{j=0}^p (-j)^{2}a_j+2\sum_{j=-1}^p (-j)^{1}b_j=1$,
• au moins d'ordre 3 il faut qu'elle soit au moins d'ordre 2 et que $\displaystyle\sum_{j=0}^p (-j)^{3}a_j+3\sum_{j=-1}^p (-j)^{2}b_j=1$,
• au moins d'ordre 4 il faut qu'elle soit au moins d'ordre 2 et que $\displaystyle\sum_{j=0}^p (-j)^{4}a_j+4\sum_{j=-1}^p (-j)^{3}b_j=1$.

Dans notre cas

• la méthode est consistante;
• pour que la méthode sois d'ordre au moins 2 il faut de plus $1= a_1+4a_2-2(b_1+2b_2-b_{-1})$ ce qui est satisfait;
• pour que la méthode sois d'ordre au moins 3 il faut de plus $1= -a_1-8a_2+3(b_1+4b_2+b_{-1})$ ce qui est satisfait;
• pour que la méthode sois d'ordre au moins 4 il faut de plus $1= a_1+16a_2-4(b_1+8b_2-b_{-1})$ ce qui n'est pas satisfait.

La méthode est donc d'ordre 3 (cf. validation avec sympy ci-dessous).

# j = 0...p
# q = p+1 # nb de pas
# ω # ordre de la méthode

CAS_GENERAL = True
p = 2

q = p+1 
ordre_max = q+2 if q%2==0 else q+1

print(f"{'★'*70}\nC'est une méthode à q = {q} pas d'ordre ω ≤ {ordre_max}")


if CAS_GENERAL:
    aa = sympy.symbols(f'a_0:{q}')
    bb = sympy.symbols(f'b_0:{q}')
    bm1 = sympy.Symbol('b_{-1}')
else : # cas particulier
    aa = [sympy.Rational(18,11), -sympy.Rational(9,11),sympy.Rational(2,11)]
    bb = [0,0,0]
    bm1 = sympy.Rational(6,11)

i=1
sa=sum( [-j*aa[j] for j in range(len(aa))] )
sb=bm1+sum( [bb[j] for j in range(len(bb))] )
print(f"{'★'*70}\nLa méthode est d'ordre ω ≥ {i} (= consistante) si ")
prlat(sum(aa).factor(),"=1" )
prlat((sa).factor()+(i*sb).factor(),"=1")

for i in range(2,ordre_max+1):
    sa=sum( [(-j)**i*aa[j] for j in range(len(aa))] )
    sb=bm1+sum( [(-j)**(i-1)*bb[j] for j in range(1,len(bb))] )
    print(f"{'★'*70}\nLa méthode est d'ordre ω ≥ {i} si elle est d'ordre {i-1} et, de plus, ")
    prlat((sa).factor()+(i*sb).factor(),"=1" )

★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★
C'est une méthode à q = 3 pas d'ordre ω ≤ 4
★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★
La méthode est d'ordre ω ≥ 1 (= consistante) si 

$\displaystyle a_{0} + a_{1} + a_{2}=1$

$\displaystyle - a_{1} - 2 a_{2} + b_{0} + b_{1} + b_{2} + b_{-1}=1$

★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★
La méthode est d'ordre ω ≥ 2 si elle est d'ordre 1 et, de plus, 

$\displaystyle a_{1} + 4 a_{2} - 2 \left(b_{1} + 2 b_{2} - b_{-1}\right)=1$

★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★
La méthode est d'ordre ω ≥ 3 si elle est d'ordre 2 et, de plus, 

$\displaystyle - a_{1} - 8 a_{2} + 3 \left(b_{1} + 4 b_{2} + b_{-1}\right)=1$

★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★
La méthode est d'ordre ω ≥ 4 si elle est d'ordre 3 et, de plus, 

$\displaystyle a_{1} + 16 a_{2} - 4 \left(b_{1} + 8 b_{2} - b_{-1}\right)=1$

Q3 [3 points] Vérifier empiriquement l'ordre de convergence sur le problème de Cauchy $$\begin{cases} y'(t) = -\dfrac{y(t)}{2}\left(1-\dfrac{y(t)}{3}\right), &\forall t \in I=[0,10],\\ y(0) = 1, \end{cases} $$

# variables globales
t0     = 0
tfinal = 10
y0     = 1

phi = lambda t,y: -y/2*(1-y/3) #-y-t**2

##############################################
# solution exacte
##############################################
t   = sympy.Symbol('t')
y   = sympy.Function('y')
edo = sympy.Eq( sympy.diff(y(t),t) , phi(t,y(t)) )
solgenLIST = sympy.dsolve(edo)
display(solgenLIST)
solgen=solgenLIST[0]
display(solgen)
consts = sympy.solve( sympy.Eq( y0, solgen.rhs.subs(t,t0)) , dict=True)[0]
display(consts)
solpar = solgen.subs(consts).simplify()
display(solpar)

sol_exacte = sympy.lambdify(t,solpar.rhs,'numpy')

##############################################
# schéma (initialisation avec sol exacte pour ordre de convergence)
##############################################
def multipas(phi, tt, sol_exacte):
    h  = tt[1] - tt[0]
    uu = [sol_exacte(tt[i]) for i in range(3)]
    for i in range(2,len(tt) - 1):
        eq = lambda x : -x+18/11*uu[i]-9/11*uu[i-1]+2/11*uu[i-2]+h*6/11*phi(tt[i+1],x)
        temp = fsolve( eq ,uu[i])
        uu.append( temp[0] )
    return uu

##############################################
# ordre
##############################################
H = []
err = []
N = 10

figure(figsize=(16,5))
ax1 = subplot(1,2,1)

for k in range(6):
    N += 20
    tt = linspace(t0, tfinal, N + 1)
    H.append( tt[1] - tt[0] )
    yy = sol_exacte(tt)
    uu = multipas(phi, tt, sol_exacte)
    err.append( norm(uu-yy,inf) )
    ax1.plot(tt,uu,label=f'Approchée avec N={N}')

ax1.plot(tt,yy,label='Exacte')
xlabel('$t$')
ylabel('$y$')
ax1.grid(True)
ax1.legend()

ax2 = subplot(1,2,2)
ax2.loglog( H, err, 'r-o')
xlabel('$h$')
ylabel('$e$')
title(f'Multipas ordre = {polyfit(log(H),log(err),1)[0]:1.2f}')
ax2.grid(True)

$\displaystyle \left[ y{\left(t \right)} = \frac{3 \left(e^{C_{1}} - \sqrt{e^{C_{1} + t}}\right)}{e^{C_{1}} - e^{t}}, \ y{\left(t \right)} = \frac{3 \left(e^{C_{1}} + \sqrt{e^{C_{1} + t}}\right)}{e^{C_{1}} - e^{t}}\right]$

$\displaystyle y{\left(t \right)} = \frac{3 \left(e^{C_{1}} - \sqrt{e^{C_{1} + t}}\right)}{e^{C_{1}} - e^{t}}$

$\displaystyle \left\{ C_{1} : - \log{\left(4 \right)}\right\}$

$\displaystyle y{\left(t \right)} = \frac{3 \left(2 \sqrt{e^{t}} - 1\right)}{4 e^{t} - 1}$