2023 Contrôle Continu – M64 – Data Science & Scientific Computing 6 : Schémas numériques pour les problèmes de Cauchy

2023 Contrôle Continu

Auteur·rice

Gloria Faccanoni

Date de publication

24 mars 2024

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

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

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

Code 
%reset -f
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# variables globales
t0     = 0
tfinal = 10

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

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

plt.figure(figsize=(10,5))

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

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

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

  1. Étudier la zéro-stabilité de cette méthode.
  2. La méthode est elle consistente? Est elle convergente? Si oui, étudier théoriquement l’ordre de convergence.
  3. 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} \)

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

Code 
%reset -f
import sympy 
 
aa_eval = [0,1]
bb_eval = [0,1]
bm1_eval = 1 


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

from IPython.display import display, Math, Markdown

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)+(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}')

display(Markdown(f"C'est une méthode à $q={p + 1}$ pas d'ordre $\\omega\\le{omega}$ avec"))

texte = ''.join([ f"{sympy.latex(aa[j])}={sympy.latex(aa_eval[j])},"+r"\quad " for j in range(p+1)])
texte += ''.join([ f"{sympy.latex(bb[j])}={sympy.latex(bb_eval[j])},"+r"\quad " for j in range(p+1)])
texte += f"{sympy.latex(bm1)}={sympy.latex(bm1_eval)}"
display(Math( texte ))  

        
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"\\xi({i}) = {sympy.latex(xi(i,aa,bb,bm1))} = {sympy.latex(xi_eval[i])}" ))
    if xi_eval[i]!=1:
        display(Markdown( f"La méthode est donc d'ordre $\\omega = {i-1}$" ))   
        break

C’est une méthode à \(q=2\) pas d’ordre \(\omega\le4\) avec

\(\displaystyle a_{0}=0,\quad a_{1}=1,\quad b_{0}=0,\quad b_{1}=1,\quad b_{-1}=1\)

\(\displaystyle \xi(0) = a_{0} + a_{1} = 1.0\)

\(\displaystyle \xi(1) = - a_{1} + b_{0} + b_{1} + b_{-1} = 1\)

\(\displaystyle \xi(2) = a_{1} - 2 b_{1} + 2 b_{-1} = 1\)

\(\displaystyle \xi(3) = - a_{1} + 3 b_{1} + 3 b_{-1} = 5\)

La méthode est donc d’ordre \(\omega = 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} \)

Code 
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 12})
from scipy.optimize import fsolve


# variables globales
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 = np.zeros_like(tt)
    uu[0:2] = sol_exacte(tt[0: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[i+1] = temp[0]
    return uu

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

plt.figure(figsize=(16,5))

ax1 = plt.subplot(1,2,1)
plt.xlabel('$t$')
plt.ylabel('$y$')
ax1.grid(True)

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 = multipas(phi, tt, sol_exacte)
    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') # sur la grille la plus fine
ax1.legend()


ax2 = plt.subplot(1,2,2)
ax2.plot( np.log(H), np.log(err), 'r-o')
plt.xlabel(r'$\ln(h)$')
plt.ylabel(r'$\ln(e)$')
plt.title(f'Multipas ordre = {np.polyfit(np.log(H),np.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} \)

  1. Déterminer les coefficients \(\sigma\), \(\vartheta\), \(\gamma\) pour que le schéma soit d’ordre 3. Est-il d’ordre 4 ?
  2. É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} \)
  3. Étudier théoriquement la A-stabilité du schéma optimal.
Code 
%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 [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\).

Code 
sigma = sympy.Symbol(r"\sigma")
theta = sympy.Symbol(r"\vartheta")
gamma = sympy.Symbol(r"\gamma")


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**"))
    matrice_Butcher(s, A, b, c)
    
    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**"))
    ordre_2(s, A, b, c, j)

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

    return None


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 None


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))
    return None


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))
    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))
    return None


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))
    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))
    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))
    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))
    return None


# APPLICATION A NOTRE CAS
c = [sigma,sympy.S(1)]
b = [gamma,1-gamma]
A = [[sigma,sympy.S(0)],[theta,1-theta]]
s = len(c)
display(Markdown(f"La méthode de Runge-Kutta est à {s} étages"))
ordre_RK(s,A,b,c)

La méthode de Runge-Kutta est à 2 étages

Matrice de Butcher

\(\displaystyle \left[\begin{matrix}\sigma & \sigma & 0\\1 & \vartheta & 1 - \vartheta\\ & \gamma & 1 - \gamma\end{matrix}\right]\)

On a 3 conditions pour avoir consistance = pour être d’ordre 1

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

\(\displaystyle \sum_{j=1}^s a_{0j}=\sigma\text{ doit être égale à }\sigma\)

\(\displaystyle \sum_{j=1}^s a_{1j}=1\text{ doit être égale à }1\)

On doit ajouter 1 condition 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}\)

On doit ajouter 2 conditions 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(\gamma - 1\right) + \left(\gamma - 1\right) \left(\vartheta - 1\right)\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=\gamma \sigma^{3} - \gamma + 1\text{ doit être égale à }\frac{1}{4}\)

\(\displaystyle \sum_{i,j=1}^s b_i c_i a_{ij} c_j=\gamma \sigma^{3} - \sigma \vartheta \left(\gamma - 1\right) + \left(\gamma - 1\right) \left(\vartheta - 1\right)\text{ doit être égale à }\frac{1}{8}\)

\(\displaystyle \sum_{i,j=1}^s b_i a_{ij} c_j^2=\gamma \sigma^{3} - \sigma^{2} \vartheta \left(\gamma - 1\right) + \left(\gamma - 1\right) \left(\vartheta - 1\right)\text{ doit être égale à }\frac{1}{12}\)

\(\displaystyle \sum_{i,j,k=1}^s b_i a_{ij} a_{jk} c_k=\gamma \sigma^{3} - \sigma^{2} \vartheta \left(\gamma - 1\right) + \sigma \vartheta \left(\gamma - 1\right) \left(\vartheta - 1\right) + \left(1 - \gamma\right) \left(\vartheta - 1\right)^{2}\text{ doit être égale à }\frac{1}{24}\)

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

Code 
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(Math(sympy.latex(cnst)))


c_substituted = [expr.subs(cnst) for expr in c]
b_substituted = [expr.subs(cnst) for expr in b]
A_substituted = [[expr.subs(cnst) for expr in row] for row in A]

\(\displaystyle \left\{ \gamma : \frac{3}{4}, \ \sigma : \frac{1}{3}, \ \vartheta : 1\right\}\)

Remarque : le choix optimal conduit au même schéma du DM de 2020, l’étude a donc déjà été fait en TD et en TP !

Le schéma n’est pas d’ordre 4 :

Code 
ordre_RK(s,A_substituted,b_substituted,c_substituted)

Matrice de Butcher

\(\displaystyle \left[\begin{matrix}\frac{1}{3} & \frac{1}{3} & 0\\1 & 1 & 0\\ & \frac{3}{4} & \frac{1}{4}\end{matrix}\right]\)

On a 3 conditions pour avoir consistance = pour être d’ordre 1

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

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

\(\displaystyle \sum_{j=1}^s a_{1j}=1\text{ doit être égale à }1\)

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

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

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

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

\(\displaystyle x < 6\)