None 2023-CC-corrige

2023 Contrôle Continu

⏱ 2h30, deposer le notebook dans Moodle

EXERCICE 1 [4 points]: choix d'un schéma

Soit le problème de Cauchy $$ \begin{cases} y'(t)=\left(1-y(t)-\dfrac{1}{4+t}\right)\sin(2\pi y(t)), &t\in[0;10]\\ y(0)=y_0 \ \in\ ]0.7;0.9[ \end{cases} $$ On ne peut pas calculer analytiquement la solution exacte.

Q [4 points] Pour calculer des solutions approchées on doit utiliser un schéma numérique. Quelle caractéristique doit-il avoir ? Justifier la réponse.

🐸 Il s'agit d'un des problèmes de Cauchy proposés lors du CT 2023 de l'ECUE M61 Équations différentielles.

Bien qu'on ne puisse pas calculer analytiquement la solution exacte, nous pouvons tracer les lignes de courant pour avoir une idée du comportement de la solution en fonction de la donnée initiale. On en déduit que le problème de Cauchy est numériquement mal posé. Il faut alors utiliser un schéma d'ordre elevé.

Pour s'en convaincre nous pouvons aussi demander à odeint du module scipy de tracer la solution approchée pour différentes valeurs de $y_0$.

In [36]:
%reset -f
%matplotlib inline

from matplotlib.pylab import *
from scipy.integrate import odeint

# variables globales
t0     = 0
tfinal = 10

phi = lambda t,y: sin(pi*2*y)*(1-1/(4+t)-y)

N = 200
tt = linspace(t0, tfinal, N + 1)

figure(figsize=(10,5))

# On trace d'abord les lignes de courant
g1  = linspace(0,10,21)
g2  = linspace(0.5,1,21)
T,Y = meshgrid(g1,g2)  
DT, DY = 1, phi(T,Y)          # compute growth rate on the grid
M = hypot(DT,DY) # = sqrt(DT**2+DY**2)         # norm growth rate 
streamplot(T,Y, DT/M, DY/M, color=M, linewidth=0.5, cmap=cm.viridis, density=2, arrowstyle='->', arrowsize=1.5)
grid()
xlabel('t')
ylabel('y')
title('Lignes de courant');

# On trace quelque solution approchée
for y0 in linspace(0.7,0.9,7):
    uu = odeint(phi,y0,tt, tfirst=True)
    plot(tt,uu,lw=4,label=f'$y_0={y0:g}$')

legend(bbox_to_anchor=(1.04, 1), loc='upper left');

Exercice 2 [7 points] : étude d'un schéma multipas

On notera $\varphi_k\stackrel{\text{déf}}{=}\varphi(t_k,u_k)$. Considérons la méthode multipas $$ u_{n+1}=u_{n-1}+h\left( \varphi_{n+1}+\varphi_{n-1}\right) $$

On remarque qu'il s'agit de la méthode de Crank-Nicolson sur l'intervalle $[t_{n-1},t_{n+1}]$.

Q2 [1 points] Étudier la zéro-stabilité.

Le schéma est de la forme $$ u_{n+1} = a_0u_n+a_1u_{n-1}+hb_0\varphi_n+hb_1\varphi_{n-1}+hb_{-1}\varphi_{n+1}. $$ On a

  • $p=1$, i.e. c'est un méthode à $q=p+1=2$ pas;
  • $a_0=0$ et $a_1=1$,
  • $b_0=0$ et $b_1=1$,
  • $b_{-1}=1\neq0$ ($\implies$ la méthode est implicite).

Le premier polynôme caractéristique est $$ \varrho(r)= %r^{p+1}-\sum_{j=0}^p a_jr^{p-j}= r^2-a_0r-a_1= r^2-1 =(r-1)(r+1). $$

Les racines sont $$ r_0=1,\quad r_{1}=-1. $$ La méthode est donc zéro-stable car $$ \begin{cases} |r_j|\le1 \quad\text{pour tout }j=0,\dots,p \\ \text{ multiplicité de $r_j$ > }1 \implies |r_j|<1 \end{cases} $$

Q2 [3 points] La méthode est elle consistente? Est elle convergente? Si oui, étudier théoriquement l'ordre de convergence.

La méthode étant zéro-stable, si elle est consistante, alors elle est convergente. Pour étudier la convergence introduisons la quantité $$ \begin{align*} \xi(i) &= \sum_{j=0}^p (-j)^{i}a_j+i\sum_{j=-1}^p (-j)^{i-1}b_j \\ &= (-0)^{i}a_0+(-1)^ia_1+i \Big( (1)^{i-1}b_{-1}+(-0)^{i-1}b_0+(-1)^{i-1}b_1 \Big) \end{align*} $$ avec $(-0)^0=1$.

Pour que la méthode soit consistante il faut que $\xi(0)=\xi(1)=1$. De plus, pour que la méthode soit d'ordre $\omega$, il faut que $\xi(i)=1$ pour tout $i\le\omega$.

De plus, $\omega$ vérifie la première barrière de Dahlquist: le schéma étant implicite à $q=2$ (pair) pas consistante et zéro-stable, $\omega=q+2\le4$.

En calculant $\xi$ on trouve que la méthode est d'ordre 2 (cf. calculs avec sympy ci-dessous).

In [37]:
%reset -f
import sympy 
 
aa_eval = [0,1]
bb_eval = [0,1]
bm1_eval = 1 

###########################################################################

from IPython.display import display, Math

def xi(i,aa,bb,bm1):
    sa = sum( [ (-j)**i * aa[j] for j in range(q) ] )
    sb = bm1+sum( [(-j)**(i-1)*bb[j] for j in range(1,q)] )
    if i==1:
        sb += bb[0]
#     return (sa).factor()+(i*sb).factor()
    return (sa)+(i*sb)

q = len(aa_eval) # nb de pas
p = q-1  # j = 0...p # coeffs du schema
omega = q+2 if q%2==0 else q+1

aa = sympy.symbols(f'a_0:{q}')
bb = sympy.symbols(f'b_0:{q}')
bm1 = sympy.Symbol('b_{-1}')

expimp = "explicite" if bm1==0 else "implicite"
print(f"{'★'*70}\nC'est une méthode {expimp} à q = {q} pas d'ordre ω ≤ {omega} avec")
texte = ""
for j in range(p+1):
        texte +=f"{sympy.latex(aa[j])}={sympy.latex(aa_eval[j])}\qquad "
texte += f"{sympy.latex(bm1)}={sympy.latex(bm1_eval)}\qquad "
for j in range(p+1):
        texte += f"{sympy.latex(bb[j])}={sympy.latex(bb_eval[j])}\qquad " 
display(Math( texte )) 
print(f"{'★'*70}")
        
xi_eval = [xi(i,aa_eval,bb_eval,bm1_eval) for i in range(omega+1)]
for i in range(omega+1):
    display(Math(f"ξ({i}) = {sympy.latex(xi(i,aa,bb,bm1))} = {sympy.latex(xi_eval[i])}") )
    if xi_eval[i]!=1:
        print(f"La méthode est d'ordre ω = {i-1}\n{'★'*70}")
        break

★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★
C'est une méthode implicite à q = 2 pas d'ordre ω ≤ 4 avec
$\displaystyle a_{0}=0\qquad a_{1}=1\qquad b_{-1}=1\qquad b_{0}=0\qquad b_{1}=1\qquad $ 
★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★
$\displaystyle ξ(0) = a_{0} + a_{1} = 1.0$
$\displaystyle ξ(1) = - a_{1} + b_{0} + b_{1} + b_{-1} = 1$
$\displaystyle ξ(2) = a_{1} - 2 b_{1} + 2 b_{-1} = 1$
$\displaystyle ξ(3) = - a_{1} + 3 b_{1} + 3 b_{-1} = 5$ 
La méthode est d'ordre ω = 2
★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★

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

In [38]:
%matplotlib inline

from matplotlib.pylab import *
rcParams.update({'font.size': 12})
from scipy.optimize import fsolve

t0     = 0
tfinal = 2
y0     = -1

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

##############################################
# solution exacte
##############################################
t   = sympy.Symbol('t')
y   = sympy.Function('y')
edo = sympy.Eq( sympy.diff(y(t),t) , phi(t,y(t)) )
solgen = sympy.dsolve(edo)
# display(solgen)
consts = sympy.solve( sympy.Eq( y0, solgen.rhs.subs(t,t0)) , dict=True)[0]
# display(consts)
solpar = solgen.subs(consts)
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(2)]
    for i in range(1,len(tt) - 1):
        eq = lambda x : -x+uu[i-1]+h*(phi(tt[i+1],x)+phi(tt[i-1],uu[i-1]))
        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.plot( log(H), log(err), 'r-o')
xlabel(r'$\ln(h)$')
ylabel(r'$\ln(e)$')
title(f'Multipas ordre = {polyfit(log(H),log(err),1)[0]:1.2f}')
ax2.grid(True)
$\displaystyle y{\left(t \right)} = \frac{- 2 t + e^{t^{2}} - 2}{t + 1}$

Exercice 3 [10 points] : étude d'un schéma RK

Considérons le schéma dont la matrice de Butcher est donnée ci-dessous: $$ \begin{array}{c|cc} \sigma & \sigma & 0 \\ 1 & \vartheta & 1-\vartheta\\ \hline & \gamma & 1-\gamma \end{array} $$
In [39]:
%reset -f
%matplotlib inline

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

import sympy 
sympy.init_printing()

Q1 [5 points] Pour quelles valeurs des paramètres le schéma est d'ordre 3? Et d'ordre 4?

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

C'est un schéma implicite à $2$ étages ($s=2$) donc $\omega\le2s=4$.

In [40]:
sigma = sympy.Symbol(r"\sigma")
theta = sympy.Symbol(r"\vartheta")
gamma = sympy.Symbol(r"\gamma")

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

c = [sigma,1]
b = [gamma,1-gamma]
A = [[sigma,0],[theta,1-theta]]
s = len(c)

# Affichage matrice de Butcher
print(f'{"Matrice de Butcher":_^80}')
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'{r"Conditions pour être d ordre 1" :_^80}')
texte = "\sum_{j=1}^s b_j =" + f"{sum(b)}"
texte += r"\text{ doit être égale à }1"
display(Math(texte))
for i in range(s):
    texte = f'i={i},\quad' +r'\sum_{j=1}^s a_{ij}=' + sympy.latex(sum(A[i]))
    texte += r"\text{ doit être égale à }c_i="+sympy.latex(c[i])
    display(Math( texte ))

print(f'{r"Conditions à ajouter pour être d ordre 2" :_^80}')
texte = '\sum_{j=1}^s b_j c_j=' +sympy.latex(sum([b[i]*c[i] for i in range(s)]))
texte += r"\text{ doit être égale à }\frac{1}{2}"
display(Math(texte))

print(f'{r"Conditions à ajouter pour être d ordre 3" :_^80}')
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))
texte = r'\sum_{i,j=1}^s b_i a_{ij} c_j='
texte = texte + sympy.latex(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]))
texte += r"\text{ doit être égale à }\frac{1}{6}"
display(Math(texte))

_______________________________Matrice de Butcher_______________________________
$\displaystyle \left[\begin{matrix}\sigma & \sigma & 0\\1 & \vartheta & 1 - \vartheta\\0 & \gamma & 1 - \gamma\end{matrix}\right]$ 
_________________________Conditions pour être d ordre 1_________________________
$\displaystyle \sum_{j=1}^s b_j =1\text{ doit être égale à }1$
$\displaystyle i=0,\quad\sum_{j=1}^s a_{ij}=\sigma\text{ doit être égale à }c_i=\sigma$
$\displaystyle i=1,\quad\sum_{j=1}^s a_{ij}=1\text{ doit être égale à }c_i=1$ 
____________________Conditions à ajouter pour être d ordre 2____________________
$\displaystyle \sum_{j=1}^s b_j c_j=\gamma \sigma - \gamma + 1\text{ doit être égale à }\frac{1}{2}$ 
____________________Conditions à ajouter pour être d ordre 3____________________
$\displaystyle \sum_{j=1}^s b_j c_j^2=\gamma \sigma^{2} - \gamma + 1\text{ doit être égale à }\frac{1}{3}$
$\displaystyle \sum_{i,j=1}^s b_i a_{ij} c_j=\gamma \sigma^{2} + \sigma \vartheta \left(1 - \gamma\right) + \left(1 - \gamma\right) \left(1 - \vartheta\right)\text{ doit être égale à }\frac{1}{6}$

On commence à resoudre le système dont les inconnues sont les paramètres $\sigma, \gamma, \vartheta$ pour avoir au moins l'ordre 3:

In [41]:
cnst = sympy.solve( [ sum([b[i]*c[i] for i in range(s)])-1/sympy.S(2),
                    sum([b[i]*c[i]**2 for i in range(s)])-1/sympy.S(3),
                    sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])-1/sympy.S(6)
        ] )[0]

display(cnst)

print(f'{"La matrice de Butcher devient":_^80}')
But = But.subs(cnst)
prlat(sympy.Matrix(But))
$\displaystyle \left\{ \gamma : \frac{3}{4}, \ \sigma : \frac{1}{3}, \ \vartheta : 1\right\}$ 
_________________________La matrice de Butcher devient__________________________
$\displaystyle \left[\begin{matrix}\frac{1}{3} & \frac{1}{3} & 0\\1 & 1 & 0\\0 & \frac{3}{4} & \frac{1}{4}\end{matrix}\right]$

Le schéma n'est pas d'ordre 4 car:

In [42]:
print(f'{r"Conditions à ajouter pour être d ordre 4" :_^80}')

texte = '\sum_{j=1}^s b_j c_j^3='
texte = texte + sympy.latex( sum([b[i]*c[i]**3 for i in range(s)]).subs(cnst) )
texte += r"\text{ doit être égale à }\frac{1}{4}"
display(Math(texte))

texte = '\sum_{i,j=1}^s b_i c_i a_{ij} c_j='
texte = texte + sympy.latex( sum([b[i]*c[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).subs(cnst) )
texte += r"\text{ doit être égale à }\frac{1}{8}"
display(Math(texte))

texte = '\sum_{i,j=1}^s b_i a_{ij} c_j^2=' 
texte = texte + sympy.latex( sum([b[i]*A[i][j]*c[j]**2 for i in range(s) for j in range(s)]).subs(cnst) )
texte += r"\text{ doit être égale à }\frac{1}{12}"
display(Math(texte))

texte = '\sum_{i,j,k=1}^s b_i a_{ij}a_{jk} c_k=' 
texte = texte + sympy.latex( 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)]).subs(cnst) )
texte += r"\text{ doit être égale à }\frac{1}{24}"
display(Math(texte))

____________________Conditions à ajouter pour être d ordre 4____________________
$\displaystyle \sum_{j=1}^s b_j c_j^3=\frac{5}{18}\text{ doit être égale à }\frac{1}{4}$
$\displaystyle \sum_{i,j=1}^s b_i c_i a_{ij} c_j=\frac{1}{9}\text{ doit être égale à }\frac{1}{8}$
$\displaystyle \sum_{i,j=1}^s b_i a_{ij} c_j^2=\frac{1}{18}\text{ doit être égale à }\frac{1}{12}$
$\displaystyle \sum_{i,j,k=1}^s b_i a_{ij}a_{jk} c_k=\frac{1}{18}\text{ doit être égale à }\frac{1}{24}$

Q2 [3 points] Étudier empiriquement l'ordre de convergence sur le problème de Cauchy $$ \begin{cases} y'(t)=\dfrac{t^2+y^2(t)}{ty(t)}, & t\in[1;5]\\ y(1)=1 \end{cases} $$

In [43]:
# variables globales
t0     = 1
tfinal = 5
y0     = 1

phi = lambda t,y: (t**2+y**2)/(t*y)

##############################################
# solution exacte
##############################################
import sympy 
sympy.init_printing()

t   = sympy.Symbol('t')
y   = sympy.Function('y')
edo = sympy.Eq( sympy.diff(y(t),t) , phi(t,y(t)) )
solgenLIST = sympy.dsolve(edo,y(t))
display(solgenLIST)
for solgen in solgenLIST: 
    #display(solgen)
    consts = sympy.solve( sympy.Eq( y0, solgen.rhs.subs(t,t0)) , dict=True)
    #display(consts)
    if consts!=[]:
        solpar = solgen.subs(consts[0]).simplify()
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)


##############################################
# solution approchée
##############################################
sol_exacte = sympy.lambdify(t,solpar.rhs,'numpy')

# Je note sigma_e etc pour sigma evaluée comme un float
sigma_e = 1/3
gamma_e = 1/(2*(1-sigma_e)) # pour ordre 2
theta_e = 1  

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*sigma_e , uu[i]+h*X[0]*sigma_e  ) , 
                           -X[1]+phi( tt[i]+h*theta_e , uu[i]+h*(X[0]*theta_e+X[1]*(1-theta_e))  )
                           ]
        K_start = phi(tt[i],uu[i])
        K1, K2 = fsolve( sys , [K_start,K_start ] )
        uu.append( uu[i]+h*(gamma_e*K1+(1-gamma_e)*K2) )
    return uu


##############################################
# ordre
##############################################
H = []
err = []

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

N = 50
for k in range(6):
    N += 50
    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.plot( log(H), log(err), 'r-o')
xlabel('$\ln(h)$')
ylabel('$\ln(e)$')
title(f'RK : ordre = {polyfit(log(H),log(err),1)[0]:1.2f}')
ax2.grid(True)
$\displaystyle \left[ y{\left(t \right)} = - t \sqrt{C_{1} + 2 \log{\left(t \right)}}, \ y{\left(t \right)} = t \sqrt{C_{1} + 2 \log{\left(t \right)}}\right]$
$\displaystyle y{\left(t \right)} = t \sqrt{C_{1} + 2 \log{\left(t \right)}}$
$\displaystyle \left\{ C_{1} : 1\right\}$
$\displaystyle y{\left(t \right)} = t \sqrt{2 \log{\left(t \right)} + 1}$

Q3 [2 points] Étudier théoriquement la A-stabilité du schéma.

La méthode appliquée à l'EDO $y'(t)=-\beta y(t)$ avec $\beta>0$ s'écrit: $$\begin{cases} u_0 = y_0 \\ K_1 = -\beta u_n-(\beta h) \sigma K_1 \leadsto K_1=\frac{-\beta}{1+(\beta h) \sigma}u_n\\ K_2 = -\beta\left(u_n+h K_1 \vartheta + h K_2 (1-\vartheta)\right)\leadsto K_2=\frac{1+\frac{(\beta h)\vartheta}{1+(\beta h) \sigma}}{1+(\beta h)(1-\vartheta)}(-\beta u_n)\\ u_{n+1} = u_n + \frac{h}{4}\left(3K_1+K_2\right) & n=0,1,\dots N-1 \end{cases} $$

En particulier, dans le cas d'ordre 3, on a $$\begin{cases} u_0 = y_0 \\ K_1 = -\beta u_n-\beta \frac{h}{3}K_1 \leadsto K_1=\frac{-3\beta}{3+\beta h}u_n\\ K_2 = -\beta\left(u_n+hK_1\right)\leadsto K_2=\frac{2\beta^2h-3\beta}{3+\beta h}u_n\\ u_{n+1} = u_n + \frac{h}{4}\left(3K_1+K_2\right) & n=0,1,\dots N-1 \end{cases} $$ L'étude ci-dessous montre que, dans ce cas, le schéma est A-stable ssi $h<\frac{6}{\beta}$.

In [44]:
sympy.var('u_n, h, β, K1, K2,')  

g = lambda y : -β*y
eq1 = sympy.Eq( K1 , g(u_n+h*K1*sigma))
eq2 = sympy.Eq( K2 , g(u_n+h*(K1*theta+K2*(1-theta))) )
sol = sympy.solve([eq1,eq2],[K1,K2])
sol = { key:value.factor() for key,value in sol.items() }
display(sol)

RHS = u_n+h*(gamma*K1+(1-gamma)*K2).factor()
RHS = RHS.subs(sol).simplify().factor()
texte = 'u_{n+1}=' + sympy.latex(RHS)
display(Math(texte)) 

print(r"Notons x=hβ>0. On trouve donc")
x = sympy.Symbol('x',positive=True)
RHS = RHS.subs(h*β,x).simplify()
texte = 'u_{n+1}=' + sympy.latex(RHS)
display(Math(texte)) 

print("Considérons le schéma d'ordre 3:")

display(cnst)

RHS = RHS.subs(cnst)
texte = r"\text{La méthode est A-stable ssi }|q(x)|<1\text{ avec }q(x)=\frac{u_{n+1}}{u_n}="
texte += sympy.latex(RHS/u_n)
display(Math(texte)) 
sympy.solve([-1<RHS/u_n,RHS/u_n<1])
$\displaystyle \left\{ K_{1} : - \frac{u_{n} β}{\sigma h β + 1}, \ K_{2} : \frac{u_{n} β \left(\sigma h β - \vartheta h β + 1\right)}{\left(\sigma h β + 1\right) \left(\vartheta h β - h β - 1\right)}\right\}$
$\displaystyle u_{n+1}=- \frac{u_{n} \left(\gamma \sigma h^{2} β^{2} - \gamma h^{2} β^{2} - \sigma \vartheta h^{2} β^{2} + \sigma h β + \vartheta h^{2} β^{2} - \vartheta h β + 1\right)}{\left(\sigma h β + 1\right) \left(\vartheta h β - h β - 1\right)}$ 
Notons x=hβ>0. On trouve donc
$\displaystyle u_{n+1}=\frac{u_{n} \left(\gamma \sigma x^{2} - \gamma x^{2} - \sigma \vartheta x^{2} + \sigma x + \vartheta x^{2} - \vartheta x + 1\right)}{\left(\sigma x + 1\right) \left(- \vartheta x + x + 1\right)}$ 
Considérons le schéma d'ordre 3:
$\displaystyle \left\{ \gamma : \frac{3}{4}, \ \sigma : \frac{1}{3}, \ \vartheta : 1\right\}$
$\displaystyle \text{La méthode est A-stable ssi }|q(x)|<1\text{ avec }q(x)=\frac{u_{n+1}}{u_n}=\frac{\frac{x^{2}}{6} - \frac{2 x}{3} + 1}{\frac{x}{3} + 1}$
Out[44]:
$\displaystyle x < 6$
In [ ]: