2024 Contrôle Terminal session 2

Auteur·rice

Gloria Faccanoni

Date de publication

27 juin 2024

⏱ 2h30, deposer le notebook dans Moodle

EXERCICE 1 : Étude de l’ ❤️ entre Roméo et Juliette

Considérons une histoire d’amour entre Roméo et Juliette. Nous définissons $R(t)$ comme l’amour de Roméo pour Juliette à l’instant $t$ et $J(t)$ comme l’amour de Juliette pour Roméo au même instant. Leur amour est modélisé par le système d’équations différentielles suivant :

$ \begin{cases} R'(t) = aR(t) + bJ(t) + f(t)\\ J'(t) = cR(t) + dJ(t) + g(t) \end{cases} $

Les paramètres $a$, $b$, $c$ et $d$ sont des constantes réelles qui décrivent les styles romantiques de Roméo et Juliette. Plus précisément, $a$ et $b$ caractérisent le style romantique de Roméo, tandis que $c$ et $d$ décrivent celui de Juliette.

Nous introduisons également une influence externe ou une fonction de perturbation $f(t)$ pour Roméo et $g(t)$ pour Juliette, pour tenir compte des influences ou perturbations externes.

Nous allons maintenant examiner comment leur dynamique amoureuse évolue dans différentes situations, avec différents types de fonctions de perturbation $f(t)$ et $g(t)$. Pour tous les tests, nous prenons $a = -1$, $b = 2$, $c = -2$, $d = 1$ et la valeur initiale est $(0.1; 0.1)$. Nous considérons l’intervalle de temps d’un mois (c’est-à-dire 30 jours).

Sans Perturbation Externe : Nous commençons par exclure les termes de perturbation externes dans le système : $f(t)=g(t)=0$ pour tout $t$. Nous examinerons le comportement périodique de leur amour.
Surprise pour Juliette : Supposons que Juliette reçoit une surprise de la part de Roméo, comme un cadeau inattendu. Cette situation peut être décrite par $g(t)$, une fonction de stimulation, qui n’a de valeur qu’à un jour donné : $ g(t) = \begin{cases} 12 & \text{ si }10<t<11\\ 0 & \text{ sinon} \end{cases} $ Nous examinerons le comportement de leur amour dans cette situation avec $f(t)=0$.
Nouvel Emploi pour Roméo : Supposons que Roméo obtient un nouvel emploi satisfaisant et est confiant dans leur vie future. Ici, $f(t)=-0.03t$ et $g(t)=0$ pour tout $t$.
Perturbation Périodique Mineure : Enfin, nous considérons une situation de perturbation périodique mineure, qui représente l’influence de toutes sortes de petits faits quotidiens. Ici, $f(t)=\sin(0.2\pi t)$ et $g(t)=0$ pour tout $t$.

EXERCICE : utiliser le schéma de Heun pour approcher la solution du système d’équations différentielles. Pour chaque situation, tracer sur le même plan les graphiques de $R(t)$ et $J(t)$ sur l’intervalle $[0, 30]$. Tracer ensuite le portrait de phase de leur amour, c’est-à-dire le graphe de $J(t)$ en fonction de $R(t)$.

Les schémas

Code

%reset -f

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

t0, R0, J0 = 0, 0.1, 0.1
tfin = 30 # on étudie l'évolution sur 30 jours

a, b, c, d = -1, 2, -2, 1
pphi = lambda t, yy, f, g : np.array([a*yy[0]+b*yy[1]+f(t), c*yy[0]+d*yy[1]+g(t)]) 

def HEUN(pphi, f, g, yy0, tt):
    RR_JJ = np.zeros((len(tt), len(yy0)))
    RR_JJ[0] = yy0
    h = tt[1]-tt[0]
    for i in range(len(tt)-1):
        k1 = pphi(tt[i], RR_JJ[i], f, g)
        k2 = pphi(tt[i+1], RR_JJ[i]+h*k1, f, g)
        RR_JJ[i+1] = RR_JJ[i] + h/2*(k1+k2)
    return RR_JJ

def EM(pphi, f, g, yy0, tt):
    RR_JJ = np.zeros((len(tt), len(yy0)))
    RR_JJ[0] = yy0
    h = tt[1]-tt[0]
    for i in range(len(tt)-1):
        k1 = pphi(tt[i], RR_JJ[i], f, g)
        RR_JJ[i+1] = RR_JJ[i] + h*pphi(tt[i]+h/2, RR_JJ[i]+h/2*k1, f, g)
    return RR_JJ    

def deja(pphi, f, g, yy0, tt):
    RR_JJ = solve_ivp(pphi, (t0, tfin), yy0, t_eval=tt, args=(f, g), 
                    #   method='BDF'
                    #   method='RK45'
                    #   method='RK23'
                    #   method='DOP853'
                    #   method='Radau'
                      method='LSODA'
                      )
    return RR_JJ.y.T

La fonction pour l’affichage

Code

def affichage(heures=24):
    mask = range(heures)
    fig = plt.figure(figsize=(24, 8))
    for i, (RR, JJ, title) in enumerate(zip([RR_test1, RR_test2, RR_test3, RR_test4], 
                                            [JJ_test1, JJ_test2, JJ_test3, JJ_test4], 
                                            ["1. Sans Perturbation Externe", "2. Surprise pour Juliette", "3. Nouvel emploi pour Roméo", "4. Perturbation périodique mineure"])):
        RR_local, JJ_local, tt_local = RR[mask], JJ[mask], tt[mask]
        plt.subplot(2, 4, i+1)
        plt.title(title)
        plt.plot(tt_local, RR_local, label='R(t)')
        plt.plot(tt_local, JJ_local, label='J(t)')
        plt.xlim([t0, tfin])
        plt.xlabel('t')
        plt.legend()
        plt.subplot(2, 4, i+1+4)
        plt.plot(RR_local, JJ_local)
        plt.scatter(RR[0], JJ[0], c='r', label='t=0')
        plt.scatter(RR_local[-1], JJ_local[-1], c='g', label=f't={tt_local[-1] : .2f}', marker='x')
        plt.xlabel('R(t)')
        plt.ylabel('J(t)')# Ajouter les axes x=0 et y=0
        plt.axhline(0, color='black',linewidth=0.5)
        plt.axvline(0, color='black',linewidth=0.5)
        plt.xlim([min(RR), max(RR)])
        plt.ylim([min(JJ), max(JJ)])
        plt.legend()

    fig.tight_layout()

Les calculs

Code

tt = np.linspace(t0, tfin, 24*tfin+1) # on veut une solution toutes les heures

schema = 'HEUN'
# schema = 'EM'
# schema = 'deja'


# 1. Sans Perturbation Externe
f = lambda t: 0
g = lambda t: 0
RR_test1, JJ_test1 = eval(schema)(pphi, f, g, [R0, J0], tt).T

# 2. Surprise pour Juliette au jour 10
f = lambda t: 0
g = lambda t: 12 if 10<t<11 else 0
RR_test2, JJ_test2 = eval(schema)(pphi, f, g, [R0, J0], tt).T

# 3. Nouvel emploi pour Roméo
f = lambda t: -0.03*t
g = lambda t: 0
RR_test3, JJ_test3 = eval(schema)(pphi, f, g, [R0, J0], tt).T

# 4. Perturbation riodique mineure
f = lambda t: np.sin(0.2*np.pi*t)
g = lambda t: 0
RR_test4, JJ_test4 = eval(schema)(pphi, f, g, [R0, J0], tt).T

Affichage statique

Code

affichage(24*11+1) # 11 jours
affichage(24*tfin+1) # 30 jours

Sans Perturbation Externe : le plan de phase suggère que leur amour est périodique. Leurs sentiments l’un pour l’autre oscillent périodiquement.
Surprise pour Juliette : le plan de phase suggère que ce type de perturbation ou de stimulus romantique peut entraîner un amour mutuel.
Nouvel Emploi pour Roméo : le plan de phase suggère que leur relation se développe vers un amour mutuel.
Perturbation Périodique Mineure : il est difficile de prévoir le développement de leur relation.

Affichage dynamique (si `interact` installé)

Code

from ipywidgets import interact
interact(affichage, heures=(24,24*tfin,24)); # on increment par 24h

EXERCICE 2 : schéma 👣 explicite avec paramètres

Pour quelles valeurs de $\eta$ et $\vartheta$ le schéma multipas linéaire suivant est convergent ? $ u_{n+1} = (3\vartheta-2\eta)u_n + \left(2-\frac{5}{2}\eta\right)u_{n-1} + h \eta(\varphi_n+\varphi_{n-1}) $
Sans effectuer de calculs, quel est l’ordre maximal théorique de ce schéma ?
Calculer l’ordre de convergence de ce schéma.
Tester empiriquement l’ordre de convergence du schéma ainsi obtenu sur le problème de Cauchy $y'(t) = \sin(2t)-y(t)\sin(t)$ et $y(\pi/2) = 1$ pour $t\in[\pi/2,2\pi]$.

Remarques :

les deux premiers points peuvent être résolus “à la main” ou avec sympy (dans ma correction, je vous montrerai comment faire avec sympy).
dans le dernier point, attention à bien indiquer à chaque étape si les fonctions que vous utlisez appartiennent à numpy ou à sympy.

Correction point 1

C’est un schéma à $q=2$ pas, c’est-à-dire de la forme $ u_{n+1} = a_0u_n+a_1u_{n-1}+h\big(b_{-1}\varphi_{n+1}+b_0\varphi_n+b_1\varphi_{n-1}\big) $ avec $a_0=3\vartheta-2\eta$, $a_1=2-\frac{5}{2}\eta$, $b_{-1}=0$ et $b_0=b_1=\eta$. Il est donc explicite.

Pour qu’il soit convergent, il faut que le schéma soit zéro-stable et consistant.

Écrivons donc qui sont $\xi(0)$, $\xi(1)$, $\xi(2)$, $\xi(3)$ puis le premier polynôme caractéristique $\varrho$ et calculons ses racines dans le cas général d’un schéma à $q=2$ pas :

Code

%reset -f

from sympy import symbols, Symbol, latex, factor, solve, Eq, poly, Rational, roots
from IPython.display import display, Math, Markdown, Latex

# ========================
q = 2
# ========================


a_0, a_1 = symbols('a_0 a_1 ', real=True)
b_0, b_1, b_m1 = symbols('b_0 b_1 b_{-1}', real=True)
r = Symbol('r')

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

# ========================
aa = [a_0, a_1]
bb = [b_0, b_1]

display(Markdown(f"Pour une méthode à q={q} pas"))   

display(Markdown(f"- on écrit le premier polynôme caractéristique (pour étudier la zéro-stabilité)"))   
rho = poly(r**(q) - sum( aa[j] * r**(q-1-j) for j in range(q)), r)
texte =  r"$\varrho(r)=" + latex(rho.as_expr(),order='old') +  r"$"
display(Math(texte))
display(Markdown(f"dont les racines sont :"))
racines = roots(rho,r)
for k,v in racines.items():
    display(k)

display(Markdown(f"- on affiche $\\xi(i)$ (pour étudier la consistance et l'ordre de consistance)"))
for i in range(4):
    display(Markdown(f"$\\xi({i}) = {latex(xi(i, q))}$"))

Pour une méthode à q=2 pas

on écrit le premier polynôme caractéristique (pour étudier la zéro-stabilité)

$\displaystyle \varrho(r)=- a_{1} - r a_{0} + r^{2}$

dont les racines sont :

$\displaystyle \frac{a_{0}}{2} - \frac{\sqrt{a_{0}^{2} + 4 a_{1}}}{2}$

$\displaystyle \frac{a_{0}}{2} + \frac{\sqrt{a_{0}^{2} + 4 a_{1}}}{2}$

on affiche $\xi(i)$ (pour étudier la consistance et l’ordre de consistance)

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

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

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

$\xi(3) = - a_{1} + 3 \left(b_{1} + b_{-1}\right)$

Étude de la consistance

Pour que le schéma soit consistant il faut que $\xi(0)=\xi(1)=1$: $ \begin{cases} a_0+a_1=1 \\ -a_1+b_0+b_1+b_{-1}=1 \end{cases} $ Dans notre cas, cela donne

Code

eta, theta = symbols(r'\eta \vartheta', real=True)

aa = [3*theta-2*eta, 2-5*eta/2]
bb = [eta, eta]
b_m1 = 0

systeme = [Eq(xi(i, q),1) for i in range(2)]
display(Markdown(r"$\begin{cases}"+' '.join([  r"\displaystyle"+latex(s)+r"\\" for s in systeme]) + r"\end{cases}$"))

$\begin{cases}\displaystyle- \frac{9 \eta - 6 \vartheta - 4}{2} = 1\\ \displaystyle2 \eta + \frac{5 \eta - 4}{2} = 1\\\end{cases}$

dont la solution est

Code

sol = solve(systeme)
for cle, valeur in sol.items():
    display(Math( latex(cle) + "=" + latex(valeur)))

$\displaystyle \eta=\frac{2}{3}$

$\displaystyle \vartheta=\frac{2}{3}$

ainsi

Code

for i in range(2):
    display(Markdown(f"$a_{i} = {latex(aa[i])} = {latex(aa[i].subs(sol))}$"))
    display(Markdown(f"$b_{i} = {latex(bb[i])} = {latex(bb[i].subs(sol))}$"))
display(Markdown(f"$b_{{-1}} = {latex(b_m1)} $"))

$a_0 = - 2 \eta + 3 \vartheta = \frac{2}{3}$

$b_0 = \eta = \frac{2}{3}$

$a_1 = 2 - \frac{5 \eta}{2} = \frac{1}{3}$

$b_1 = \eta = \frac{2}{3}$

$b_{-1} = 0 $

Étude de la zéro-stabilité

Pour que le schéma soit zéro-stable il faut que le premier polynôme caractéristique $\varrho(r)=r^2-a_0r-a_1$ ait ses racines de module inférieur ou égal à 1. Si une racine est de module 1, elle doit être simple.

Le premier polynôme caractéristique est $r^2+2r-3$. Ses racines sont $r_1=-1/3$ et $r_2=1$. Le schéma est donc zéro-stable.

Code

rho = poly(r**(q) - sum( aa[j].subs(sol) * r**(q-1-j) for j in range(q)), r)
texte =  r"$\varrho(r)=" + latex(rho.as_expr(),order='old') + r"$"
display(Math(texte))
racines = roots(rho.subs(sol),r)
display(racines)

$\displaystyle \varrho(r)=- \frac{1}{3} - \frac{2 r}{3} + r^{2}$

{1: 1, -1/3: 1}

Correction point 2

D’après la première barrière de Dahlquist, un schéma à $q=2$ pas explicite est au mieux d’ordre $q=2$. On doit donc juste vérifier si $\xi(2)=1$.

Code

xi2 = xi(2, q)
texte = r"\xi(2) =" 
texte += latex(xi2)
texte += "=" + latex(xi2.subs(sol))

display(Math(texte))

$\displaystyle \xi(2) =- 2 \eta - \frac{5 \eta - 4}{2}=-1$

On conclut que le schéma n’est que d’ordre 1.

Correction point 3

Dans cette partie, attention à bien indiquer à chaque étape si les fonctions utlisées appartiennent à numpy ou à sympy.

Code

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


##################################################
# EDO sympy
##################################################
y0     = 1
t0_sym = sympy.pi/2
phi_sym = lambda t,z: -sympy.sin(t)*z+sympy.sin(2*t)

##################################################
# solution exacte sympy
##################################################
t, z   = sympy.symbols('t z')
y   = sympy.Function('y')
edo = sympy.Eq( sympy.diff(y(t),t) , phi_sym(t,y(t)) )
solpar = sympy.dsolve(edo, ics={y(t0_sym): y0})
display(solpar)

##################################################
# sympy -> numpy 
##################################################

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

t0, tfinal = float(t0_sym), 2*np.pi
y0 = float(y0)

phi = sympy.lambdify((t,z),phi_sym(t,z),'numpy')

##################################################
# schéma 
##################################################
def multipas(phi, tt, sol_exacte):
    h  = tt[1] - tt[0]
    
    uu = np.zeros_like(tt)
    
    q = 2
    a0, a1,  = 2/3, 1/3 
    b0, b1 = 2/3, 2/3

    # initialisation avec solution exacte pour ordre de convergence
    uu[0:q+1] = sol_exacte(tt[0:q+1])
    # recurrence    
    for i in range(q, len(tt) - 1):
        uu[i+1] = a0*uu[i] + a1*uu[i-1] + h*(b0*phi(tt[i],uu[i]) + b1*phi(tt[i-1],uu[i-1]) )
    return uu

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


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

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

N = 100
for k in range(6):
    N += 10
    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'Le schéma a ordre {np.polyfit(np.log(H),np.log(err),1)[0]:1.2f}')
ax2.grid(True)

$\displaystyle y{\left(t \right)} = - e^{\cos{\left(t \right)}} + 2 \cos{\left(t \right)} + 2$

EXERCICE 3 : schéma 🎨 RK à 3 étages explicite

La matrice de Butcher d’un schéma RK à $s=3$ étages quelconque s’écrit : $ \begin{array}{c|ccc} c_0 & a_{00} & a_{01} & a_{02} \\ c_1 & a_{10} & a_{11} & a_{12} \\ c_2 & a_{20} & a_{21} & a_{22} \\ \hline & b_0 & b_1 & b_2 \end{array} $

Écrire la matrice de Butcher correspondante au schéma suivant : $ \begin{cases} u_0 = y_0 \\ k_0 = \varphi\Big(t_n,u_n\Big) \\ k_1 = \varphi\Big(t_n+\frac{h}{3},u_n+\frac{h}{3}k_0\Big) \\ k_2 = \varphi\Big(t_n+\frac{2h}{3},u_n+\frac{2h}{3}k_1\Big) \\ u_{n+1}= u_n + h\left(b_0k_0+(1-b_0)k_2\right) \end{cases} $
Déterminer pour quelle valeur de $b_0$ l’ordre de convergence de ce schéma est maximal.
Tester empiriquement l’ordre de convergence de ce schéma sur le problème de Cauchy $y'(t) = 4t(t^2-1)+y(t)$ et $y(0) = 1$ pour $t\in[0,1]$.
Déterminer sous quelle condition sur $h$ ce schéma est A-stable.

Correction point 1 : matrice de Butcher

$ \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\(1em] \frac{1}{3} & \frac{1}{3} & 0 & 0 \\(1em] \frac{2}{3} & 0 & \frac{2}{3} & 0 \\(1em] \hline & b_0 & 0 & (1-b_0) \end{array} $

Il s’agit d’un schéma explicite, car la matrice $\mathbb A$ est triangulaire inférieure stricte.

Correction point 2 : borne supérieure pour l’ordre et calcul de $b_0$

La barrière de Dalquist affirme qu’un schéma RK à $s=3$ étages explicite ne peut pas être d’ordre supérieur à $s=3$.

On peut bien-sûr écrire les formules à la main, mais on peut aussi utiliser un solveur pour résoudre le système d’équations non-linéaires comme suit :

Code

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


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**"))
    Eq_1 = ordre_1(s, A, b, c, j)
    Eqs.append(Eq_1)
    
    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.append(Eq_3)  
    
    # display(Markdown("**On doit ajouter 4 conditions pour être d'ordre 4**"))
    # Eq_4 = ordre_4(s, A, b, c, j)
    # Eqs.append(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))
    Eq_1 = [sympy.Eq( sum(b).simplify(), 1 )]
    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 ))
        Eq_1.append(sympy.Eq( somma, c[i] ))
    return Eq_1    


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

s = 3

# paramètres dans le cas générale
b_0, b_1, b_2 = sympy.symbols('b_0 b_1 b_2', real=True)
c_0, c_1, c_2 = sympy.symbols('c_0 c_1 c_2', real=True)
a_00, a_01, a_02, a_10, a_11, a_12, a_20, a_21, a_22 = sympy.symbols('a_00 a_01 a_02 a_10 a_11 a_12 a_20 a_21 a_22', real=True)

# matrice dans le cas général
b = [b_0, b_1, b_2]
c = [c_0, c_1, c_2]
A = [[a_00, a_01, a_02], [a_10, a_11, a_12], [a_20, a_21, a_22]]

# matrice dans notre cas particulier
b = [ b_0, 0, 1-b_0]
c = [ 0, sympy.Rational(1,3), sympy.Rational(2,3) ]
A = [[ 0, 0, 0], [ sympy.Rational(1,3), 0, 0], [ 0, sympy.Rational(2,3), 0]]

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

Matrice de Butcher

$\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\\frac{1}{3} & \frac{1}{3} & 0 & 0\\\frac{2}{3} & 0 & \frac{2}{3} & 0\\ & b_{0} & 0 & 1 - b_{0}\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}=\frac{2}{3}\text{ doit être égale à }\frac{2}{3}$

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

$\displaystyle \sum_{j=1}^s b_j c_j=\frac{2}{3} - \frac{2 b_{0}}{3}\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{4}{9} - \frac{4 b_{0}}{9}\text{ doit être égale à }\frac{1}{3}$

$\displaystyle \sum_{i,j=1}^s b_i a_{ij} c_j=\frac{2}{9} - \frac{2 b_{0}}{9}\text{ doit être égale à }\frac{1}{6}$

En conclusion, le schéma est d’ordre 1 par construction. De plis, le paramètre $b_0$ doit vérifier les conditions suivantes pour que le schéma soit d’ordre 2 et éventuellement 3 :

Code

for i,eq in enumerate(Eqs[1:]):
    display(Markdown(f"**Pour être d'ordre {i+2} il faut vérifier**"))
    display(Math(sympy.latex((eq))))

Pour être d’ordre 2 il faut vérifier

$\displaystyle \left[ \frac{2}{3} - \frac{2 b_{0}}{3} = \frac{1}{2}\right]$

Pour être d’ordre 3 il faut vérifier

$\displaystyle \left[ \frac{4}{9} - \frac{4 b_{0}}{9} = \frac{1}{3}, \ \frac{2}{9} - \frac{2 b_{0}}{9} = \frac{1}{6}\right]$

La première condition qui garantie l’ordre 2 impose le choix de $b_0$ :

Code

sol = sympy.solve(Eqs[1])   
sol

$\displaystyle \left\{ b_{0} : \frac{1}{4}\right\}$

Vérifions si cette valeur satisfait les deux conditions pour être d’ordre 3:

Code

Eqs[2][0].subs(sol), Eqs[2][1].subs(sol)

$\displaystyle \left( \text{True}, \ \text{True}\right)$

On conclut que le schéma avec $b_0=\frac{1}{4}$ est d’ordre 3.

Code

# Le vecteur b est donc
b = [ b_0.subs(sol), 0, (1-b_0).subs(sol),]

Correction point 3 : test de l’ordre

Code

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



##################################################
# EDO
##################################################
t0     = 0
tfinal = 1
y0     = 1

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

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

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

##################################################
# schéma 
##################################################

# Paramètres du schéma recuperés depuis la matrice de Butcher obtenue précédemment
cc = np.array(c)
bb = np.array(b)
AA = np.array(A)

# Schéma de Runge-Kutta explicite à 4 étages
def RK(phi, tt, y0):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    # initialisation
    uu[0] = y0
    # recurrence    
    for i in range(len(tt) - 1):
        k0 = phi(tt[i]          , uu[i] )
        k1 = phi(tt[i] + cc[1]*h, uu[i] + h* AA[1,0]*k0)
        k2 = phi(tt[i] + cc[2]*h, uu[i] + h*(AA[2,0]*k0 + AA[2,1]*k1))
        uu[i+1] = uu[i] + h*(bb[0]*k0 + bb[1]*k1 + bb[2]*k2)
    return uu


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


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

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

N = 10
for k in range(6):
    N += 10
    tt = np.linspace(t0, tfinal, N + 1)
    H.append( tt[1] - tt[0] )
    yy = sol_exacte(tt)
    uu = RK(phi, tt, 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') # 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'Le schéma a ordre {np.polyfit(np.log(H),np.log(err),1)[0]:1.2f}')
ax2.grid(True)

$\displaystyle y{\left(t \right)} = - 4 t^{3} - 12 t^{2} - 20 t + 21 e^{t} - 20$

Correction point 4 : A-stabilité

On applique le schéma au problème de Cauchy $y'(t) = -\beta y(t)$ et $y(0) = 1$ avec $\beta>0$ dont la solution est $y(t)=e^{-\beta t}$.

On obtient alors la relation de récurrence suivante :

Code

u_n = sympy.symbols(r'u_n')
u_np1 = sympy.symbols(r'u_{n+1}')
beta = sympy.symbols(r'\beta', positive=True)
dt = sympy.symbols(r'h', positive=True)
x = sympy.symbols(r'x', positive=True)

k0 = -beta*u_n
k1 = -beta*(u_n + dt* AA[1,0]*k0)
k2 = -beta*(u_n + dt*(AA[2,0]*k0 + AA[2,1]*k1))
u_np1 = u_n + dt*(bb[0]*k0 + bb[1]*k1 + bb[2]*k2)

display(Math(sympy.latex(sympy.Eq(sympy.Symbol('u_{n+1}'),u_np1.simplify()))))

$\displaystyle u_{n+1} = \frac{u_{n} \left(- \beta h \left(\beta h \left(\beta h - 3\right) + 6\right) + 6\right)}{6}$

On peut expliciter $u_n$ en fonction de $u_0$ et du produit $\beta h$ (on notera $x$ le produit $\beta h$) :

Code

q = (u_np1/u_n).simplify().subs(beta*dt,x)
display(Math(sympy.latex(sympy.Eq(sympy.Symbol('q(x)'),q.expand().simplify()))))

$\displaystyle q(x) = - \frac{x^{3}}{6} + \frac{x^{2}}{2} - x + 1$

La suite vérifie $u_n \xrightarrow[n\to \infty]{} 0$ ssi $\left|q(x)\right|<1$. En étudiant brièvement la fonction $q(x)$, on trouve que $q(x)<1$ ssi $x<x_0$ donné ci-dessous:.

Code

sympy.plot(q, (x,0,5), ylim=(-1.5,1.5),
            ylabel='q(x)', xlabel='x', title='q(x) en fonction de x')

sol = sympy.solveset(q+1, x, domain=sympy.S.Reals)  
display(Markdown(f"$\\beta h<{sol.args[0].evalf()}$"))

$\beta h<2.51274532661833$

EXERCICE 1 : Étude de l’ ❤️ entre Roméo et Juliette

Les schémas

La fonction pour l’affichage

Les calculs

Affichage statique

Affichage dynamique (si interact installé)

EXERCICE 2 : schéma 👣 explicite avec paramètres

Correction point 1

Correction point 2

Correction point 3

EXERCICE 3 : schéma 🎨 RK à 3 étages explicite

Correction point 1 : matrice de Butcher

Correction point 2 : borne supérieure pour l’ordre et calcul de \(b_0\)

Correction point 3 : test de l’ordre

Correction point 4 : A-stabilité

Affichage dynamique (si `interact` installé)