None 2023-CT-corrige

2023 Contrôle Terminal session 1

⏱ 2h30, deposer le notebook dans Moodle

EXERCICE 1 [7 points] : étude d'un système

On s'intèresse au problème de Cauchy $$ \begin{cases} y'(t)=\sin(y(t))-2z(t)e^{y(t)},\\ z'(t)=z(t)\Big(z(t)e^{y(t)}-\cos(y(t))\Big),\\ y(0)=\frac{1}{2},\\ z(0)=\frac{1}{4}, \end{cases} \qquad t\in[0,30] $$

Q1 [1 point] Notons $$E(y,z) = z\sin(y)-z^2e^y$$ et $\mathcal{E}(t)=E(y(t),z(t))$. Montrer que $\mathcal{E}$ est un invariant.

On remarque que $\frac{\partial E}{\partial y}(y(t),z(t))=-z'(t)$ et $\frac{\partial E}{\partial z}(y(t),z(t))=y'(t)$. On a alors $$ \mathcal{E}'(t)= \frac{d }{dt}E(y(t),z(t))= \frac{\partial E}{\partial y}y'(t)+\frac{\partial E}{\partial z}z'(t)= -z'(t)y'(t)+y'(t)z'(t)=0. $$ Donc $\mathcal{E}(t)=\mathcal{E}(0)=0$.

Q2 [2 points] Calculer la solution approchée obtenue avec le schéma d'Euler Explicite avec $N=100$. Dans trois repères différents afficher

  • dans le premier repère $t\mapsto y(t)$ et $t\mapsto z(t)$ dans le même,
  • dans le deuxième repère $y\mapsto z(y)$,
  • dans le troisième repère $t\mapsto\mathcal{E}(t)$.
In [66]:
%reset -f
%matplotlib inline

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

def affichage(tt,uu,ww,HH_num,s):

    suptitle(f'Schéma {s}')
    #tight_layout()

    subplot(1,3,1)
    plot(tt,uu,'o--',tt,ww,'d--')
    xlabel('t')
    legend([r'$y(t)$',r'$z(t)$'])
    grid()

    subplot(1,3,2)
    plot(uu,ww,'o--')
    xlabel('y')
    ylabel('z')
    grid()
    Y1,Y2 = meshgrid(linspace(min(uu),max(uu),21),linspace(min(ww),max(ww),21))
    
    V1,V2 = pphi(tt,[Y1,Y2]) 
    r1 = sqrt(1+V1**2)
    r2 = sqrt(1+V2**2)
    quiver(Y1, Y2, V1/r1, V2/r2, cmap=cm.viridis, scale_units='xy')
    # axis('equal')

    subplot(1,3,3)
    plot(tt,HH_num,'o--')
    xlabel('t')
    ylabel(r'$\mathcal{E}$')
    grid();

    
t0 = 0
y0, z0 = 0.5, 0.25

tfinal = 30
N = 100

phi1 = lambda t,y,z : sin(y)-2*z*exp(y)
phi2 = lambda t,y,z : z*(z*exp(y)-cos(y))
Ham = lambda y,z : z*sin(y)-z**2*exp(y) 

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

# pour odeint
pphi = lambda t,yy : [phi1(t,yy[0],yy[1]),phi2(t,yy[0],yy[1])]
yy0  = [y0,z0]
In [67]:
def EE(phi1,phi2,tt,y0,z0):
    h  = tt[1]-tt[0]
    uu = [y0]
    ww = [z0]
    for i in range(len(tt)-1):
        uu.append(uu[i]+h*phi1(tt[i],uu[i],ww[i]))
        ww.append(ww[i]+h*phi2(tt[i],uu[i],ww[i])) 
    return [uu,ww]


[uu, ww] = EE(phi1,phi2,tt,y0,z0)
HH_num = [ Ham(u,w) for (u,w) in zip(uu,ww) ]

figure(figsize=(18,5))
affichage(tt,uu,ww,HH_num,"EE")

Q3 [3 points]
Même exercice avec le schéma de Cranck Nicolson.

In [68]:
def CN(phi1,phi2,tt,y0,z0):
    h  = tt[1]-tt[0]
    uu = [y0]
    ww = [z0]
    for i in range(len(tt)-1):
        sys = lambda vv : [-vv[0]+uu[i]+h/2*(phi1(tt[i],uu[i],ww[i])+phi1(tt[i+1], vv[0], vv[1])), 
                           -vv[1]+ww[i]+h/2*(phi2(tt[i],uu[i],ww[i])+phi2(tt[i+1], vv[0], vv[1])) ]
        sol_u, sol_w = fsolve( sys, [uu[i],ww[i]] )
        uu.append(sol_u)
        ww.append(sol_w)
    return [uu,ww]

[uu, ww] = CN(phi1,phi2,tt,y0,z0)
HH_num = [ Ham(u,w) for (u,w) in zip(uu,ww) ]

figure(figsize=(18,5))
affichage(tt,uu,ww,HH_num,"CN")

Q4 [1 points]
Même exercice avec odeint.

In [69]:
from scipy.integrate import odeint

uu,ww = odeint(pphi,yy0,tt,tfirst=True).T
HH_num = [ Ham(u,w) for (u,w) in zip(uu,ww) ]

figure(figsize=(18,5))
affichage(tt,uu,ww,HH_num,"scipy.odeint")

Exercice 2 [8 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}=\frac{4}{3}u_n-\frac{1}{3}u_{n-1}+h\frac{2}{3}\varphi_{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}+h\left(b_0\varphi_n+b_1\varphi_{n-1}\right)+hb_{-1}\varphi_{n+1}. $$ On a

  • $p=1$: c'est un méthode à $q=p+1=2$ pas;
  • $a_0=\frac{4}{3}$ et $a_1=-\frac{1}{3}$,
  • $b_0=0$ et $b_1=0$,
  • $b_{-1}=\frac{2}{3}\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-\frac{4}{3}r+\frac{1}{3} =(r-1)\left(r-\frac{1}{3}\right). $$

Les racines sont $$ r_0=1,\quad r_{1}=\frac{1}{3}. $$ 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 [4 points] La méthode est elle consistente? Est elle convergente? Si oui, étudier théoriquement l'ordre de convergence.

On a prouvé que la méthode est zéro-stable. Ainsi, si elle est consistante, alors elle est convergente.

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) \\ &= (-0)^{i}\frac{4}{3}-(-1)^i\frac{1}{3}+ (1)^{i-1}i \frac{2}{3} \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:

  • $\xi(0) = (-0)^{0}\frac{4}{3}-(-1)^0\frac{1}{3}+ (1)^{0-1}0 \frac{2}{3} = \frac{4}{3}-\frac{1}{3} = 1$

  • $\xi(1) = (-0)^{1}\frac{4}{3}-(-1)^1\frac{1}{3}+ (1)^{1-1}1 \frac{2}{3} = \frac{1}{3}+ \frac{2}{3}=1$

  • $\xi(2) = (-0)^{2}\frac{4}{3}-(-1)^2\frac{1}{3}+ (1)^{2-1}2 \frac{2}{3} = -\frac{1}{3}+ 2 \frac{2}{3} = 1$

  • $\xi(3)=\frac{1}{3}+ 3 \frac{2}{3}=\frac{7}{3}$

(vérification avec sympy ci-dessous).

In [70]:
%reset -f

import sympy 
 
aa_eval = [sympy.Rational(4,3), -sympy.Rational(1,3)]
bb_eval = [0,0]
bm1_eval = sympy.Rational(2,3) 

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

from IPython.display import display, Math, Latex

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

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


print(f"{'★'*70}\nC'est une méthode à q = {q} pas d'ordre ω ≤ {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 ))     
print(f"{'★'*70}\n")
        
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 à q = 2 pas d'ordre ω ≤ 4 avec
$\displaystyle a_{0}=\frac{4}{3},\quad a_{1}=- \frac{1}{3},\quad b_{0}=0,\quad b_{1}=0,\quad b_{-1}=\frac{2}{3}$ 
★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★
$\displaystyle ξ(0) = a_{0} + a_{1} = 1$
$\displaystyle ξ(1) = - a_{1} + b_{0} + b_{1} + b_{-1} = 1$
$\displaystyle ξ(2) = a_{1} - 2 \left(b_{1} - b_{-1}\right) = 1$
$\displaystyle ξ(3) = - a_{1} + 3 \left(b_{1} + b_{-1}\right) = \frac{7}{3}$ 
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^2}, &\forall t \in I=[2,5],\\ y(2) = 1, \end{cases} $$

In [71]:
%matplotlib inline

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


# variables globales
t0     = 2
tfinal = 5
y0     = 1

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

##############################################
# 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+4/3*uu[i]-1/3*uu[i-1]+h*2/3*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.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{C_{1} e^{- 2 t}}{\left(t - 1\right)^{\frac{3}{2}} \sqrt{t + 1}} - 2$
$\displaystyle y{\left(t \right)} = -2 + \frac{3 \sqrt{3} e^{4} e^{- 2 t}}{\left(t - 1\right)^{\frac{3}{2}} \sqrt{t + 1}}$

Exercice 3 [8 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|ccc} 0 & 0 & 0 & 0 \\ \eta & \eta & 0 & 0\\ \sigma & 0 & \sigma & 0\\ \hline & 1-\gamma & 0 & \gamma \end{array} $$ **Dans le sujet de cet exercice, il y a un paramètre $\eta$. Pour le fixer, remplacez mon nom et prénom par les vôtres et vous obtiendrez la valeur pour ce paramètre dans tout l'exercice.**
In [72]:
%reset -f
%matplotlib inline
from IPython.display import display, Math
import sympy 
sympy.init_printing()

# =======================================================
nom = "Faccanoni"
prenom = "Gloria"
# =======================================================
L=list(range(2,5))
idx=sum([ord(c) for c in nom+prenom])%len(L)
print("Paramètre de l'exercice ")
display(Math(f"\\eta = {sympy.latex(sympy.Rational(1,L[idx]))}"))

Paramètre de l'exercice
$\displaystyle \eta = \frac{1}{3}$

Q1 [3 points] Pour quelles valeurs des paramètres $\sigma$ et $\gamma$ le schéma est d'ordre 2? Et d'ordre 3?

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

  • C'est un schéma explicite à $3$ étages ($s=3$) donc $\omega\le s=3$.

  • Il est consistant (i.e. d'ordre 1) pour tout $\sigma$ et tout $\gamma$.

Les conditions pour être d'ordre 2 et 3 sont rappelées ci-dessous. On trouve que :

  • il est d'ordre 2 ssi $\sigma\gamma=\frac{1}{2}$,
  • il est d'ordre 3 ssi, de plus, $\sigma^2\gamma=\frac{1}{3}$ et $\frac{\sigma\gamma\eta}{3}=\frac{1}{6}$.

Conclusion : il est au moins d'ordre 2 ssi $\sigma\gamma=\frac{1}{2}$, et d'ordre 3 ssi $\sigma=\frac{2}{3}$, $\gamma=\frac{3}{4}$ et $\eta=\frac{1}{3}$.

In [73]:
from matplotlib.pylab import *
from scipy.optimize import fsolve

import sympy 
sympy.init_printing()

sigma = sympy.Symbol(r"\sigma")
gamma = sympy.Symbol(r"\gamma")

# =======================================================
eta = sympy.Rational(1,3)
# =======================================================

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

c = [0,eta,sigma]
b = [1-gamma,0,gamma]
A = [[0,0,0],[eta,0,0],[0,sigma,0]]
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}0 & 0 & 0 & 0\\\frac{1}{3} & \frac{1}{3} & 0 & 0\\\sigma & 0 & \sigma & 0\\0 & 1 - \gamma & 0 & \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}=\mathtt{\text{0}}\text{ doit être égale à }c_i=0$
$\displaystyle i=1,\quad\sum_{j=1}^s a_{ij}=\frac{1}{3}\text{ doit être égale à }c_i=\frac{1}{3}$
$\displaystyle i=2,\quad\sum_{j=1}^s a_{ij}=\sigma\text{ doit être égale à }c_i=\sigma$ 
____________________Conditions à ajouter pour être d ordre 2____________________
$\displaystyle \sum_{j=1}^s b_j c_j=\gamma \sigma\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}\text{ doit être égale à }\frac{1}{3}$
$\displaystyle \sum_{i,j=1}^s b_i a_{ij} c_j=\frac{\gamma \sigma}{3}\text{ doit être égale à }\frac{1}{6}$

Q2 [3 points] Étudier empiriquement l'ordre de convergence sur le problème de Cauchy $$ \begin{cases} y'(t)=\left(\cos(2\pi t)-4t^2\right)y(t), & t\in[0;5]\\ y(0)=4 \end{cases}$$ On affichera les courbes de convergence pour un schéma d'ordre 1 (à vous de donner un choix de $\sigma$ et $\gamma$ pour obtenir un tel schéma) et pour un schéma d'ordre 2 (ou 3 si possible).

💭 ne pas confondre sympy.exp ou sympy.pi ... (du module sympy pour le calcul de la solution exacte) et exp pu pi ... (du module numpy pour le calcul numérique).

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

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


# variables globales
t0     = 0
tfinal = 5
y0     = 4

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

phi_sym = lambda t,y: (4*sympy.cos(2*sympy.pi*t)-4*t**2)*y

t   = sympy.Symbol('t')
y   = sympy.Function('y')
edo = sympy.Eq( sympy.diff(y(t),t) , phi_sym(t,y(t)) )
solgen = sympy.dsolve(edo,y(t))
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')

phi = lambda t,y: (4*cos(2*pi*t)-4*t**2)*y


def RK(tt,phi,y0,sigma,gamma,eta):
    uu = [y0]
    h = tt[1]-tt[0]
    for i in range(len(tt)-1):
        K1 = phi( tt[i] , uu[i] )
        K2 = phi( tt[i]+h*eta , uu[i]+h*eta*K1 )
        K3 = phi( tt[i]+h*sigma , uu[i]+h*sigma*K2 )
        uu.append( uu[i]+h*((1-gamma)*K1+gamma*K3) )
    return uu

##############################################
# ordre
##############################################
H = []
err_1, err_2, err_3 = [], [], []

figure(figsize=(16,10))
ax11 = subplot(3,2,1)
ax12 = subplot(3,2,2)
ax21 = subplot(3,2,3)
ax22 = subplot(3,2,4)
ax31 = subplot(3,2,5)
ax32 = subplot(3,2,6)

ax11.set_title("Schéma d'ordre 1")
ax21.set_title("Schéma d'ordre 2")
ax32.set_title("Schéma d'ordre 3")

N = 0
for k in range(6):
    N += 100
    tt = linspace(t0, tfinal, N + 1)
    H.append( tt[1] - tt[0] )
    yy = sol_exacte(tt)

    sigma, gamma, eta = 1/2, 1/2, 1/3 # cas ordre 1
    uu_1 = RK(tt,phi,y0,sigma,gamma,eta)
    err_1.append( norm(uu_1-yy,inf) )
    ax11.plot(tt,uu_1,label=f'Approchée avec N={N}')   
    
    sigma, gamma, eta = 1/4, 2, 1/6 # cas ordre 2
    uu_2 = RK(tt,phi,y0,sigma,gamma,eta)
    err_2.append( norm(uu_2-yy,inf) )
    ax21.plot(tt,uu_2,label=f'Approchée avec N={N}')   
    
    sigma, gamma, eta = 2/3, 3/4, 1/3 # cas ordre 3
    uu_3 = RK(tt,phi,y0,sigma,gamma,eta)
    err_3.append( norm(uu_3-yy,inf) )
    ax31.plot(tt,uu_3,label=f'Approchée avec N={N}')    
    

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

ax12.plot( log(H), log(err_1), 'r-o')
xlabel('$\ln(h)$')
ylabel('$\ln(e)$')
ax12.set_title(f'RK : ordre = {polyfit(log(H),log(err_1),1)[0]:1.2f}')
ax12.grid(True)

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

ax22.plot( log(H), log(err_2), 'r-o')
xlabel('$\ln(h)$')
ylabel('$\ln(e)$')
ax22.set_title(f'RK : ordre = {polyfit(log(H),log(err_2),1)[0]:1.2f}')
ax22.grid(True)


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

ax32.plot( log(H), log(err_3), 'r-o')
xlabel('$\ln(h)$')
ylabel('$\ln(e)$')
ax32.set_title(f'RK : ordre = {polyfit(log(H),log(err_3),1)[0]:1.2f}')
ax32.grid(True)
$\displaystyle y{\left(t \right)} = C_{1} e^{- \frac{4 t^{3}}{3} + \frac{2 \sin{\left(2 \pi t \right)}}{\pi}}$
$\displaystyle \left\{ C_{1} : 4\right\}$
$\displaystyle y{\left(t \right)} = 4 e^{- \frac{4 t^{3}}{3} + \frac{2 \sin{\left(2 \pi t \right)}}{\pi}}$

Q3 [2 points] Étudier théoriquement la A-stabilité du schéma d'ordre maximal (dans la correction ci-dessous $\eta=\frac{1}{3}$ et donc l'ordre maximal est 3).

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\\ K_2 = -\beta \left(u_n+\frac{h}{3}K_1\right)\leadsto K_2=\left(1-\frac{\beta h}{3}\right)K_1\\ K_3 = -\beta \left(u_n+\frac{2h}{3}K_2\right)\leadsto K_3=\left(1-\frac{2\beta h}{3}\left(1-\frac{\beta h}{3}\right)\right)K_1\\ u_{n+1} = u_n + \frac{h}{4}\left(K_1+4K_3\right) & n=0,1,\dots N-1 \end{cases} $$ L'étude ci-dessous montre que le schéma est A-stable ssi $h<\frac{\vartheta}{\beta}$ avec $\vartheta\approx2.7$.

In [75]:
%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('u_n, h, β, K1, K2, sigma, gamma, eta')  
g = lambda y : -β*y
K1 = g(u_n)
K2 = g(u_n+h*eta*K1)
K3 = g(u_n+h*sigma*K2)
RHS = u_n+h*((1-gamma)*K1+gamma*K3)
RHS = RHS.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)) 
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)) 

# ======== On remplace les paramètres par les valeurs optimales ==========
RHS = RHS.subs(gamma,sympy.Rational(4,3)).subs(sigma,sympy.Rational(2,3)).subs(eta,sympy.Rational(1,3)).simplify()
# ======================================================================

texte = r"\text{Dans le cas optimal on a } \sigma=\frac{2}{3}, \gamma=\frac{4}{3} \text{ ainsi } \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])


q = sympy.apart(RHS/u_n)
texte = r"\text{On pose }" + "q(x):="  + sympy.latex( q )
display( Math( texte ) )

prlat( r'q(0)=' , q.subs(x,0) )

dq=sympy.diff(q,x).simplify()
prlat( r"q'(x)=" , dq )

sol=sympy.solveset(dq,domain=sympy.S.Reals)
texte = r"q'(x)<0 \text{ pour tout } x \text{ et l'on a} "
display( Math( texte ) )

sol = sympy.solveset(q+1,domain=sympy.S.Reals)
texte = r"q(x)=-1 \iff x= " + sympy.latex(sol) + r"\approx" + sympy.latex(sol.evalf()) 
display( Math( texte ) )

from matplotlib.pylab import *
xx=linspace(0,5,100)
yy=array([q.subs(x,xi).evalf() for xi in xx])
plot(xx,yy)
plot([xx[0],xx[-1]],[-1,-1])
plot([xx[0],xx[-1]],[1,1])
axis([0,5,-2,2])
grid();
$\displaystyle u_{n+1}=- u_{n} \left(\eta \gamma h^{3} \sigma β^{3} - \gamma h^{2} \sigma β^{2} + h β - 1\right)$ 
Notons x=hβ>0. On trouve donc
$\displaystyle u_{n+1}=u_{n} \left(- \eta \gamma \sigma x^{3} + \gamma \sigma x^{2} - x + 1\right)$
$\displaystyle \text{La méthode est A-stable ssi }|q(x)|<1\text{ avec }q(x)=\frac{u_{n+1}}{u_n}=- \eta \gamma \sigma x^{3} + \gamma \sigma x^{2} - x + 1$
$\displaystyle \text{Dans le cas optimal on a } \sigma=\frac{2}{3}, \gamma=\frac{4}{3} \text{ ainsi } \frac{u_{n+1}}{u_n}=- \frac{8 x^{3}}{27} + \frac{8 x^{2}}{9} - x + 1$
$\displaystyle \text{On pose }q(x):=- \frac{8 x^{3}}{27} + \frac{8 x^{2}}{9} - x + 1$
$\displaystyle \mathtt{\text{q(0)=}}1$
$\displaystyle \mathtt{\text{q'(x)=}}- \frac{8 x^{2}}{9} + \frac{16 x}{9} - 1$
$\displaystyle q'(x)<0 \text{ pour tout } x \text{ et l'on a} $
$\displaystyle q(x)=-1 \iff x= \left\{- \frac{1}{24 \sqrt[3]{\frac{43}{432} + \frac{\sqrt{822}}{288}}} + 1 + 3 \sqrt[3]{\frac{43}{432} + \frac{\sqrt{822}}{288}}\right\}\approx\left\{2.68038081801272\right\}$