2022 Contrôle Terminal session 2

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

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

1. Calculer et tracer le graphe de la solution exacte \(t\mapsto y(t)\) d’abord pour \(b=g=10\) puis pour \(b=g=100\).

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

   • 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.
   • Soit \(b=g=100\). Tracer la solution exacte et les solutions approchées obtenues avec \(N=10\), \(N=50\), \(N=100\).

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

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve

import sympy 
sympy.init_printing()

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

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

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

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

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

Q2 [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. 1. Soit \(b=g=100\). 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 = np.zeros_like(tt)
    uu[0] = y0
    h = tt[1]-tt[0]
    for n in range(len(tt)-1):
        K1 = (np.exp(-b*h)-1)/(-b*h)
        uu[n+1] = np.exp(-b*h)*uu[n]+g*h*K1 
    return uu


# ======================================================
# CALCUL SOLUTION APPROCHEE
# ======================================================
plt. figure(1,figsize=(20,7))
for N in [10,50,100]:
    tt = np.linspace(t0,tfinal,N+1)
    uu = EEx(tt,y0)
    plt.plot(tt,uu,'-D',label=f'{N=}')    

yy = exacte(tt,b,b,1+1.e-8)
plt.plot(tt,yy,'-',label='Exacte')
plt.legend()
plt.grid()

2 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} \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}& \color{red}\vartheta & 0\\ 1 & \frac{3}{2} & -\frac{3}{2}& \frac{1}{2} & \frac{1}{2}\\ \hline & \frac{3}{2} & \color{red}\gamma & \frac{1}{2} & \frac{1}{2} \end{array} \)

1. Calculer \(\gamma\) et \(\vartheta\) pour que le schéma soit consistant. Dans la suite on posera \(\gamma\) et \(\vartheta\) égales à ces valeurs.

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

3. Étudier théoriquement l’ordre de convergence du schéma.

%reset -f
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve

import sympy 
sympy.init_printing()

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

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 \(\color{red}\gamma=-\dfrac{3}{2}\).
• La deuxième condition est satisfaite ssi $ = - ++ $, il faut donc poser $= $.

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

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

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

##############################################
# solution exacte
##############################################
t   = sympy.Symbol('t')
y   = sympy.Function('y')
edo = sympy.Eq( sympy.diff(y(t),t) , phi(t,y(t)) )
solpar = sympy.dsolve(edo, y(t), ics={y(t0):y0}).simplify()
display(solpar)

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

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


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

_, (ax1,ax2) = plt.subplots(nrows=1,ncols=2,figsize=(16,5))

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

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

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

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

def ordre_RK(s, A=None, b=None, c=None):

    j = sympy.symbols('j')
    
    if A is None:
        A = sympy.MatrixSymbol('a', s, s)
    else:
        A = sympy.Matrix(A)
    if c is None:
        c = sympy.symbols(f'c_0:{s}')    
    if b is None:
        b = sympy.symbols(f'b_0:{s}')
    
    
    display(Markdown("**Matrice de Butcher**"))
    But = matrice_Butcher(s, A, b, c)
    
    Eqs = []

    # display(Markdown(f"**On a {s+1} conditions pour avoir consistance = pour être d'ordre 1**"))
    # ordre_1(s, A, b, c, j)
    
    display(Markdown("**On doit ajouter 1 condition pour être d'ordre 2**"))
    Eq_2 = ordre_2(s, A, b, c, j)
    Eqs.append(Eq_2)

    display(Markdown("**On doit ajouter 2 conditions pour être d'ordre 3**"))
    Eq_3 = ordre_3(s, A, b, c, j)  
    Eqs.extend(Eq_3)  
    
    display(Markdown("**On doit ajouter 4 conditions pour être d'ordre 4**"))
    Eq_4 = ordre_4(s, A, b, c, j)
    Eqs.extend(Eq_4)

    return Eqs


def matrice_Butcher(s, A, b, c):
    But = sympy.Matrix(A)
    But = But.col_insert(0, sympy.Matrix(c))   
    last = [sympy.Symbol(" ")]
    last.extend(b)
    But = But.row_insert(s,sympy.Matrix(last).T)
    display(Math(sympy.latex(sympy.Matrix(But))))
    return But


def ordre_1(s, A, b, c, j):
    texte = "\sum_{j=1}^s b_j =" + f"{sum(b).simplify()}"
    texte += r"\text{ doit être égale à }1"
    display(Math(texte))
    for i in range(s):
        somma = sympy.summation(A[i,j],(j,0,s-1)).simplify()
        texte = r'\sum_{j=1}^s a_{'+str(i)+r'j}=' + sympy.latex( somma )
        texte += r"\text{ doit être égale à }"+sympy.latex(c[i])
        display(Math( texte ))
    return None    


def ordre_2(s, A, b, c, j):
    texte = '\sum_{j=1}^s b_j c_j=' +sympy.latex(sum([b[i]*c[i] for i in range(s)]).simplify())
    texte += r"\text{ doit être égale à }\frac{1}{2}"
    display(Math(texte))
    Eq_2 = sympy.Eq( sum([b[i]*c[i] for i in range(s)]).simplify(), sympy.Rational(1,2) )
    return Eq_2


def ordre_3(s, A, b, c, j):
    texte = '\sum_{j=1}^s b_j c_j^2='
    texte += sympy.latex( sum([b[i]*c[i]**2 for i in range(s)]).simplify() )
    texte += r"\text{ doit être égale à }\frac{1}{3}"
    display(Math(texte))
    Eq_3_1 = sympy.Eq( sum([b[i]*c[i]**2 for i in range(s)]).simplify(), sympy.Rational(1,3) )
    texte = r'\sum_{i,j=1}^s b_i a_{ij} c_j='    
    somma = sum([b[i]*A[i,j]*c[j] for j in range(s) for i in range(s)]).simplify()
    texte = texte + sympy.latex(somma)
    texte += r"\text{ doit être égale à }\frac{1}{6}"
    display(Math(texte))
    Eq_3_2 = sympy.Eq( sum([b[i]*A[i,j]*c[j] for j in range(s) for i in range(s)]).simplify(), sympy.Rational(1,6) )
    return [Eq_3_1, Eq_3_2]


def ordre_4(s, A, b, c, j):
    texte = '\sum_{j=1}^s b_j c_j^3='
    texte += sympy.latex( sum([b[i]*c[i]**3 for i in range(s)]).simplify() )
    texte += r"\text{ doit être égale à }\frac{1}{4}"
    display(Math(texte))
    Eq_4_1 = sympy.Eq( sum([b[i]*c[i]**3 for i in range(s)]).simplify(), sympy.Rational(1,4) )
    texte = r'\sum_{i,j=1}^s b_i c_i a_{ij} c_j='
    somma = sum([b[i]*c[i]*A[i,j]*c[j] for j in range(s) for i in range(s)]).simplify()
    texte = texte + sympy.latex(somma)
    texte += r"\text{ doit être égale à }\frac{1}{8}"
    display(Math(texte))
    Eq_4_2 = sympy.Eq( sum([b[i]*c[i]*A[i,j]*c[j] for j in range(s) for i in range(s)]).simplify(), sympy.Rational(1,8) )
    texte = r'\sum_{i,j=1}^s b_i a_{ij} c_j^2='
    somma = sum([b[i]*A[i,j]*c[j]**2 for j in range(s) for i in range(s)]).simplify()
    texte = texte + sympy.latex(somma)
    texte += r"\text{ doit être égale à }\frac{1}{12}"
    display(Math(texte))
    Eq_4_3 = sympy.Eq( sum([b[i]*A[i,j]*c[j]**2 for j in range(s) for i in range(s)]).simplify(), sympy.Rational(1,12) )
    texte = r'\sum_{i,j,k=1}^s b_i a_{ij} a_{jk} c_k='
    somma = sum([b[i]*A[i,j]*A[j,k]*c[k] for k in range(s) for j in range(s) for i in range(s)]).simplify()
    texte = texte + sympy.latex(somma)
    texte += r"\text{ doit être égale à }\frac{1}{24}"
    display(Math(texte))
    Eq_4_4 = sympy.Eq( sum([b[i]*A[i,j]*A[j,k]*c[k] for k in range(s) for j in range(s) for i in range(s)]).simplify(), sympy.Rational(1,24) )
    return [Eq_4_1, Eq_4_2, Eq_4_3, Eq_4_4]



##################################################
# Conditions avec s étages
##################################################

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)


Eqs = ordre_RK(s, A, b, c)

# for eq in Eqs:
#     display(eq)

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}\\ & \frac{3}{2} & - \frac{3}{2} & \frac{1}{2} & \frac{1}{2}\end{matrix}\right]\)

On doit ajouter 1 condition pour être d’ordre 2

\(\displaystyle \sum_{j=1}^s b_j c_j=\frac{1}{2}\text{ doit être égale à }\frac{1}{2}\)

On doit ajouter 2 conditions pour être d’ordre 3

\(\displaystyle \sum_{j=1}^s b_j c_j^2=\frac{1}{3}\text{ doit être égale à }\frac{1}{3}\)

\(\displaystyle \sum_{i,j=1}^s b_i a_{ij} c_j=\frac{1}{6}\text{ doit être égale à }\frac{1}{6}\)

On doit ajouter 4 conditions pour être d’ordre 4

\(\displaystyle \sum_{j=1}^s b_j c_j^3=\frac{11}{36}\text{ doit être égale à }\frac{1}{4}\)

\(\displaystyle \sum_{i,j=1}^s b_i c_i a_{ij} c_j=\frac{5}{48}\text{ doit être égale à }\frac{1}{8}\)

\(\displaystyle \sum_{i,j=1}^s b_i a_{ij} c_j^2=\frac{5}{72}\text{ doit être égale à }\frac{1}{12}\)

\(\displaystyle \sum_{i,j,k=1}^s b_i a_{ij} a_{jk} c_k=\frac{1}{48}\text{ doit être égale à }\frac{1}{24}\)

D’après ces résultats, le schéma est d’ordre \(3\).

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

1. Parmi ces deux méthode, indiquer laquelle est consistente (justifier la réponse). Dans la suite on n’étudiera que cette méthode.

2. Étudier la zéro-stabilité de la méthode.

3. Étudier théoriquement l’ordre de la méthode.

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

%reset -f
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve

import sympy as sp
sp.init_printing()

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

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} \xi(0)\stackrel{\text{def}}{=}\displaystyle\sum_{j=0}^p a_j=a_0+a_1=1, \\ \xi(1)\stackrel{\text{def}}{=}\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.

aa = [sp.Rational(18,11), sp.Rational(-9,11), sp.Rational(2,11)]
q = len(aa)

r = sp.Symbol('r')
rho = sp.Function(r'\varrho')
p = q-1
rho_gen =  sp.poly(r**(p+1) - sum( aa[j] * r**(p-j) for j in range(p+1)), r)
eq = sp.Eq( rho(r), rho_gen.as_expr() )

display(Markdown(f"Le premier polynome caractéristique est :"))
display(Math( sp.latex(eq,order='old') ) )

display(Markdown(f"Ses racines sont :"))
racines = sp.solve(rho_gen,r)
display(racines)
display(Markdown(f"dont les modules valent :"))
display([abs(rac).evalf() for rac in racines])

Le premier polynome caractéristique est :

\(\displaystyle \varrho{\left(r \right)} = - \frac{2}{11} + \frac{9 r}{11} - \frac{18 r^{2}}{11} + r^{3}\)

Ses racines sont :

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

dont les modules valent :

\(\displaystyle \left[ 1.0, \ 0.426401432711221, \ 0.426401432711221\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 
★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★
La méthode est d'ordre ω ≥ 2 si elle est d'ordre 1 et, de plus, 
★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★
La méthode est d'ordre ω ≥ 3 si elle est d'ordre 2 et, de plus, 
★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★
La méthode est d'ordre ω ≥ 4 si elle est d'ordre 3 et, de plus, 

\(\displaystyle a_{0} + a_{1} + a_{2}\mathtt{\text{=1}}\)

\(\displaystyle - a_{1} - 2 a_{2} + b_{0} + b_{1} + b_{2} + b_{-1}\mathtt{\text{=1}}\)

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

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

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

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