2024 Contrôle Terminal session 2

⏱ 2h30, deposer le notebook dans Moodle

import numpy as np
import matplotlib.pyplot as plt
import sympy as sp
from scipy.optimize import fsolve
from scipy.integrate import solve_ivp
from IPython.display import display, Markdown, Math

EXERCICE 1 : Étude de l’amour 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).

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

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

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

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

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

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

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

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

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

2. Surprise pour Juliette : le plan de phase suggère que ce type de perturbation ou de stimulus romantique peut entraîner un amour mutuel.

3. Nouvel Emploi pour Roméo : le plan de phase suggère que leur relation se développe vers un amour mutuel.

4. Perturbation Périodique Mineure : il est difficile de prévoir le développement de leur relation.

Affichage dynamique (si interact installé)

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

EXERCICE 2 : Schéma MPL explicite avec paramètres

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

2. Sans effectuer de calculs, quel est l’ordre maximal théorique de ce schéma ?
   Calculer l’ordre de convergence de ce schéma.

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

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 :

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


a_0, a_1 = sp.symbols('a_0 a_1 ', real=True)
b_0, b_1, b_m1 = sp.symbols('b_0 b_1 b_{-1}', real=True)
r = sp.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 = sp.poly(r**(q) - sum( aa[j] * r**(q-1-j) for j in range(q)), r)
texte =  r"$\varrho(r)=" + sp.latex(rho.as_expr(),order='old') +  r"$"
display(Math(texte))
display(Markdown(f"dont les racines sont :"))
racines = sp.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}) = {sp.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

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

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

systeme = [sp.Eq(xi(i, q),1) for i in range(2)]
display(Markdown(r"$\begin{cases}"+' '.join([  r"\displaystyle"+sp.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

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

\(\displaystyle \eta=\frac{2}{3}\)

\(\displaystyle \vartheta=\frac{2}{3}\)

ainsi

for i in range(2):
    display(Markdown(f"$a_{i} = {sp.latex(aa[i])} = {sp.latex(aa[i].subs(sol))}$"))
    display(Markdown(f"$b_{i} = {sp.latex(bb[i])} = {sp.latex(bb[i].subs(sol))}$"))
display(Markdown(f"$b_{{-1}} = {sp.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.

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

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

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

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

xi2 = xi(2, q)
texte = r"\xi(2) =" 
texte += sp.latex(xi2)
texte += "=" + sp.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.

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

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

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

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

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

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

phi = sp.lambdify((t,z),phi_sp(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} \)

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

2. Déterminer pour quelle valeur de \(b_0\) l’ordre de convergence de ce schéma est maximal.

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

4. Déterminer sous quelle condition sur \(h\) ce schéma est A-stable.

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

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 :

prlat = lambda *args : display(Math(''.join( sp.latex(arg) for arg in args )))


def ordre_RK(s, A=None, b=None, c=None):

    j = sp.symbols('j')
    
    if A is None:
        A = sp.MatrixSymbol('a', s, s)
    else:
        A = sp.Matrix(A)
    if c is None:
        c = sp.symbols(f'c_0:{s}')    
    if b is None:
        b = sp.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 = sp.Matrix(A)
    But = But.col_insert(0, sp.Matrix(c))   
    last = [sp.Symbol(" ")]
    last.extend(b)
    But = But.row_insert(s,sp.Matrix(last).T)
    display(Math(sp.latex(sp.Matrix(But))))
    return But


def ordre_1(s, A, b, c, j):
    texte = r"\sum_{j=1}^s b_j =" + f"{sum(b).simplify()}"
    texte += r"\text{ doit être égale à }1"
    display(Math(texte))
    Eq_1 = [sp.Eq( sum(b).simplify(), 1 )]
    for i in range(s):
        somma = sp.summation(A[i,j],(j,0,s-1)).simplify()
        texte = r'\sum_{j=1}^s a_{'+str(i)+r'j}=' + sp.latex( somma )
        texte += r"\text{ doit être égale à }"+sp.latex(c[i])
        display(Math( texte ))
        Eq_1.append(sp.Eq( somma, c[i] ))
    return Eq_1    


def ordre_2(s, A, b, c, j):
    texte = r'\sum_{j=1}^s b_j c_j=' +sp.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 = sp.Eq( sum([b[i]*c[i] for i in range(s)]).simplify(), sp.Rational(1,2) )
    return Eq_2


def ordre_3(s, A, b, c, j):
    texte = r'\sum_{j=1}^s b_j c_j^2='
    texte += sp.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 = sp.Eq( sum([b[i]*c[i]**2 for i in range(s)]).simplify(), sp.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 + sp.latex(somma)
    texte += r"\text{ doit être égale à }\frac{1}{6}"
    display(Math(texte))
    Eq_3_2 = sp.Eq( sum([b[i]*A[i,j]*c[j] for j in range(s) for i in range(s)]).simplify(), sp.Rational(1,6) )
    return [Eq_3_1, Eq_3_2]


def ordre_4(s, A, b, c, j):
    texte = r'\sum_{j=1}^s b_j c_j^3='
    texte += sp.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 = sp.Eq( sum([b[i]*c[i]**3 for i in range(s)]).simplify(), sp.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 + sp.latex(somma)
    texte += r"\text{ doit être égale à }\frac{1}{8}"
    display(Math(texte))
    Eq_4_2 = sp.Eq( sum([b[i]*c[i]*A[i,j]*c[j] for j in range(s) for i in range(s)]).simplify(), sp.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 + sp.latex(somma)
    texte += r"\text{ doit être égale à }\frac{1}{12}"
    display(Math(texte))
    Eq_4_3 = sp.Eq( sum([b[i]*A[i,j]*c[j]**2 for j in range(s) for i in range(s)]).simplify(), sp.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 + sp.latex(somma)
    texte += r"\text{ doit être égale à }\frac{1}{24}"
    display(Math(texte))
    Eq_4_4 = sp.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(), sp.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 = sp.symbols('b_0 b_1 b_2', real=True)
c_0, c_1, c_2 = sp.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 = sp.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, sp.Rational(1,3), sp.Rational(2,3) ]
A = [[ 0, 0, 0], [ sp.Rational(1,3), 0, 0], [ 0, sp.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 :

for i,eq in enumerate(Eqs[1:]):
    display(Markdown(f"**Pour être d'ordre {i+2} il faut vérifier**"))
    display(Math(sp.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\) :

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

{b_0: 1/4}

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

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

(True, True)

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

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

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

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

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

sol_exacte = sp.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\)

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 :

u_n = sp.symbols(r'u_n')
u_np1 = sp.symbols(r'u_{n+1}')
beta = sp.symbols(r'\beta', positive=True)
dt = sp.symbols(r'h', positive=True)
x = sp.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(sp.latex(sp.Eq(sp.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\)) :

q = (u_np1/u_n).simplify().subs(beta*dt,x)
display(Math(sp.latex(sp.Eq(sp.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:.

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

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

\(\beta h<2.51274532661833\)