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

%reset -f
%matplotlib inline

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

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

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

    plt.subplot(1,3,1)
    plt.plot(tt,uu,'o--',label=r'$y(t)$')
    plt.plot(tt,ww,'d--',label=r'$z(t)$')
    plt.xlabel('t')
    plt.legend()
    plt.grid()

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

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

    
t0 = 0
y0, z0 = 0.5, 0.25

tfinal = 30
N = 100

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

tt = np.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]

def EE(phi1,phi2,tt,y0,z0):
    h  = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    ww = np.zeros_like(tt)
    uu[0] = y0
    ww[0] = z0
    for i in range(len(tt)-1):
        uu[i+1] = uu[i]+h*phi1(tt[i],uu[i],ww[i])
        ww[i+1] = ww[i]+h*phi2(tt[i],uu[i],ww[i])
    return [uu,ww]


[uu, ww] = EE(phi1,phi2,tt,y0,z0)
HH_num = Ham(uu,ww)

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

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

def CN(phi1,phi2,tt,y0,z0):
    h  = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    ww = np.zeros_like(tt)
    uu[0] = y0
    ww[0] = 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])) ]
        uu[i+1], ww[i+1] = fsolve( sys, [uu[i],ww[i]] )
    return [uu,ww]

[uu, ww] = CN(phi1,phi2,tt,y0,z0)
HH_num = Ham(uu,ww)

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

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

from scipy.integrate import odeint

uu,ww = odeint(pphi,yy0,tt,tfirst=True).T
HH_num = Ham(uu,ww)

plt.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 [3 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).

%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, 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).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}')

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}=\frac{4}{3},\quad a_{1}=- \frac{1}{3},\quad b_{0}=0,\quad b_{1}=0,\quad b_{-1}=\frac{2}{3}\)

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

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

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

\(\displaystyle \xi(3) = - a_{1} + 3 \left(b_{1} + b_{-1}\right) = \frac{7}{3}\)

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^2}, &\forall t \in I=[2,5],\\ y(2) = 1, \end{cases} \)

import numpy as np
import matplotlib.pyplot as plt
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 = np.zeros_like(tt)
    uu[0:2] = sol_exacte(tt[0:2])
    # 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[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)
ax1.legend()

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


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)} = -2 + \frac{3 \sqrt{3} e^{4} e^{- 2 t}}{\left(t - 1\right)^{\frac{3}{2}} \sqrt{t + 1}}\)

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

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.

%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 \(\sigma\gamma\eta=\frac{1}{6}\); autrement dit, ssi \(\sigma\gamma=\frac{1}{2}\), \(\sigma=\frac{2}{3}\) et \(\eta=\frac{1}{3}\).

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

Ci-dessous la vérification avec sympy.

import numpy as np
import matplotlib.pyplot as plt
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)
# eta = sympy.Symbol(r"\eta")
# =======================================================

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


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 = [0,eta,sigma]
b = [1-gamma,0,gamma]
A = [[0,0,0],[eta,0,0],[0,sigma,0]]
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 à 3 étages

Matrice de Butcher

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

On a 4 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}=0\text{ doit être égale à }0\)

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

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

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

\(\displaystyle \sum_{j=1}^s b_j c_j=\gamma \sigma\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}\text{ doit être égale à }\frac{1}{3}\)

\(\displaystyle \sum_{i,j=1}^s b_i a_{ij} c_j=\sigma {b}_{2} {c}_{1} + \frac{{b}_{1} {c}_{0}}{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.cos ou sympy.pi … (du module sympy pour le calcul de la solution exacte) et np.cos ou np.pi … (du module numpy pour le calcul numérique).

%reset -f
%matplotlib inline


import numpy as np
import matplotlib.pyplot as plt
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*np.cos(2*np.pi*t)-4*t**2)*y


def RK(tt,phi,y0,sigma,gamma,eta):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    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[i+1] = uu[i]+h*((1-gamma)*K1+gamma*K3)
    return uu

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

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

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

N = 0
for k in range(6):
    N += 100
    tt = np.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( np.linalg.norm(uu_1-yy,np.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( np.linalg.norm(uu_2-yy,np.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( np.linalg.norm(uu_3-yy,np.inf) )
    ax31.plot(tt,uu_3,label=f'Approchée avec N={N}')    
    

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

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

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

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


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

ax32.plot( np.log(H), np.log(err_3), 'r-o')
plt.xlabel('$\ln(h)$')
plt.ylabel('$\ln(e)$')
ax32.set_title(f'RK : ordre = {np.polyfit(np.log(H),np.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\).

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


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

display(Markdown(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"La méthode est A-stable ssi $|q(x)|<1$ avec $q(x)=\dfrac{u_{n+1}}{u_n}="
texte += sympy.latex(RHS/u_n) + "$"
display(Markdown(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()
q = sympy.apart(RHS/u_n)
# ======================================================================

texte = r"Dans le cas optimal on a $\sigma=\dfrac{2}{3}$ et $\gamma=\dfrac{4}{3}$ ainsi \(q(x)=\dfrac{u_{n+1}}{u_n}="
texte += sympy.latex(RHS/u_n) + "\)"
display(Markdown(texte)) 

display(Markdown(r"Étudions la fonction $q(x)$ pour $x>0$:"))

sympy.solve( [ -1<RHS/u_n , RHS/u_n<1 ] )

display(Markdown(f"$q(0)={q.subs(x,0)}$"))

ell = sympy.Limit(q,x,sympy.oo)
display(Math(f"{sympy.latex(ell)}={sympy.latex(ell.doit())}"))

dq = sympy.diff(q,x).simplify()
display(Markdown(f"$\displaystyle q'(x)={sympy.latex(dq)}$"))

sol=sympy.solveset(dq,domain=sympy.S.Reals)
display(Markdown(f"$q'(x)<0$ pour tout $x>0$ et l'on a"))

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

# p1 = sympy.plot(q,(x,0,5),ylim=(-2,2),xlabel=r'$x$',ylabel=r'$q(x)$',show=False,grid=True)
# p2 = sympy.plot(1,(x,0,5),ylim=(-2,2),show=False,line_color='r')
# p3 = sympy.plot(-1,(x,0,5),ylim=(-2,2),show=False,line_color='r')
# p1.extend(p2)
# p1.extend(p3)
# p1.show()

plt.figure(figsize=(12,6))
q = sympy.lambdify(x,RHS/u_n,'numpy')
xx = np.linspace(0,5,100)
plt.plot(xx,q(xx),label=r'$q(x)$')
plt.hlines(1,0,5,'r')
plt.hlines(-1,0,5,'r')
plt.axis([0,5,-2,2])
plt.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)\)

La méthode est A-stable ssi \(|q(x)|<1\) avec \(q(x)=\dfrac{u_{n+1}}{u_n}=- \eta \gamma \sigma x^{3} + \gamma \sigma x^{2} - x + 1\)

Dans le cas optimal on a \(\sigma=\dfrac{2}{3}\) et \(\gamma=\dfrac{4}{3}\) ainsi \(q(x)=\dfrac{u_{n+1}}{u_n}=- \frac{8 x^{3}}{27} + \frac{8 x^{2}}{9} - x + 1\)

Étudions la fonction \(q(x)\) pour \(x>0\):

\(q(0)=1\)

\(\displaystyle \lim_{x \to \infty}\left(- \frac{8 x^{3}}{27} + \frac{8 x^{2}}{9} - x + 1\right)=-\infty\)

\(\displaystyle q'(x)=- \frac{8 x^{2}}{9} + \frac{16 x}{9} - 1\)

\(q'(x)<0\) pour tout \(x>0\) 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\}\)