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

2019 Contrôle Terminale session 2

Auteur·rice

Gloria Faccanoni

Date de publication

06 mars 2026

Code 
import numpy as np
import matplotlib.pyplot as plt
import sympy as sp

from IPython.display import display, Markdown, Math

EXERCICE 1 : Choix d’un schéma

Considérons le problème de Cauchy \(\begin{cases} y'(t)=y(t)-2e^{-t},\\ y(0)=a \end{cases}\) 1. Calculer la solution exacte de ce problème de Cauchy en fonction de \(a\) pour \(t>0\). 2. Tracer son graphe sur l’intervalle \([0;4]\) pour \(a=1\), \(a=1.1\) et \(a=0.9\). 3. Soit \(a=1\). Parmi les schémas EE (Euler Explicite), EI (Euler Implicite) et EM (Euler modifié), lequel faut-il choisir ? Justifier la réponse en 5 lignes au plus.

Correction

Il s’agit d’une EDO linéaire d’ordre 1, on peut donc calculer la solution à la main: \( a(t)y'(t)+b(t)y(t)=g(t) \) avec \(a(t)=1\), \(b(t)=-1\) et \(g(t)=-2e^{-t}\) donc + \(A(t)=\int \frac{b(t)}{a(t)}\mathrm{d}t=\int -1\mathrm{d}t= -t\), + \(B(t)=\int \frac{g(t)}{a(t)} e^{A(t)}\mathrm{d}t =\int -2e^{-t} e^{-t}\mathrm{d}t = e^{-2t}\),

d’où \(y(t)=ce^{-A(t)}+B(t)e^{-A(t)}=c e^{t}+ e^{-t}\).

En imposant la condition \(a=y(0)\) on trouve l’unique solution du problème de Cauchy donné: \( \boxed{y(t)=(a-1)e^{t}+e^{-t}.} \)

Sinon, on peut utiliser le module de calcul symbolique.

Code 
sp.var('t,a')
u = sp.Function('u')
f = u(t)-2*sp.exp(-t)
edo = sp.Eq(sp.diff(u(t),t),f)
print("EDO:")
display(edo)

solgen = sp.dsolve(edo,u(t)).simplify()
print("Solution générale:")
display(solgen)

t0 = 0
u0 = a 
consts = sp.solve( [solgen.rhs.subs(t,t0)-u0 ], dict=True)[0]
print(f"Prise en compte de la condition initiale u({t0})={u0}:")
display(consts)

solpar = solgen.subs(consts)
print("Solution du problème de Cauchy:")
display(solpar)

print("Affichage de trois solutions")
p = sp.plot(solpar.rhs.subs({a:1}), (t,0,4), line_color='b', legend = True, show=False)
p.extend(sp.plot(solpar.rhs.subs({a:1.1}), (t,0,4), line_color='m', legend = True, show=False))
p.extend(sp.plot(solpar.rhs.subs({a:0.9}), (t,0,4), line_color='r', legend = True, show=False))
p[0].label="a=1" 
p[1].label="a=1.1" 
p[2].label="a=0.9"
p.show();
EDO:

\(\displaystyle \frac{d}{d t} u{\left(t \right)} = u{\left(t \right)} - 2 e^{- t}\)

Solution générale:

\(\displaystyle u{\left(t \right)} = C_{1} e^{t} + e^{- t}\)

Prise en compte de la condition initiale u(0)=a:
{C1: a - 1}
Solution du problème de Cauchy:

\(\displaystyle u{\left(t \right)} = \left(a - 1\right) e^{t} + e^{- t}\)

Affichage de trois solutions

Lorsque \(a=1\), le problème est numériquement mal posé, i.e. sensible aux perturbations sur la donnée initiale.
Toute méthode numérique induisant une erreur sur la solution, cette erreur va s’amplifier avec les itérations succéssives. La seule issue pour traiter ce genre de problème est d’utiliser des schémas d’ordre élévé (ici donc le schéma EM qui est d’ordre 2 tandis que les schémas EE et EI sont d’ordre 1) avec un pas \(h\) petit.

EXERCICE 2 : Circuit LC

Un circuit LC est un circuit électrique contenant une bobine (L) et un condensateur (Capacité). C’est ainsi qu’on obtient le phénomène de résonance électrique.

LC_circuit.png

Dans un circuit électrique de type LC en série, le courant \(I\), en ampères (A), et la charge électrique du condensateur \(q\), en coulombs (C), évoluent avec le temps selon le système d’équations différentielles \( \begin{cases} q'(t)=I(t)\\ L I'(t)+ \dfrac{1}{C}q(t)=0\\ \end{cases} \) avec

  • \(L>0\) l’inductance de la bobine, en henrys (H) ;
  • \(C>0\) la capacité électrique du condensateur, en farads (F) ;
  • \(t\) le temps, en secondes (s).

Pour nos calculs, on prendra les valeurs suivantes :

  • \(L=2\) H,
  • \(C=0.5\) F,
  • \(q(0)=1\) C,
  • \(i(0)=2\) A.
  1. Écrire le système sous forme matricielle ainsi que sous la forme d’une EDO d’ordre 2 en l’inconnue \(t\mapsto q(t)\). Calculer la solution exacte et afficher \(t\mapsto q(t)\), \(t\mapsto i(t)\) et \(q\mapsto i(q)\).
  2. Afficher \(t\mapsto q(t)\), \(t\mapsto i(t)\) et \(q\mapsto i(q)\) en comparant solution exacte et solution approchée obtenue avec la méthode d’Euler Explicite pour \(t\in[0;7]\) avec \(N=101\) points.
  3. Même exercice pour la méthode de Heun.
  4. Même exercice pour la méthode de Euler-Cromer. Cette méthode, appelée Euler-Cromer ou Euler symplectique ou Newton-Störmer-Verlet, s’écrit comme suit : \(\begin{cases} I_{n+1}=I_n+h\varphi_2\big(t_{n},q_{n},I_{n}\big)\\ q_{n+1}=q_{n}+h\varphi_1(t_{n+1},q_{n+1},I_{n+1}), \qquad\leftarrow\text{ce calcul est explicite car $\varphi_1$ ne dépend pas de $q$} \end{cases}\)

Correction Q1

On a l’écriture matricielle \( \begin{pmatrix} q \\ I \end{pmatrix}'(t) {}= \begin{pmatrix} 0 & 1\\ -\frac{1}{LC}&0 \end{pmatrix} \begin{pmatrix} q\\ I \end{pmatrix}(t) \) soit encore \( \begin{cases} q'(t)=\varphi_1(t,q(t),I(t)) \\ I'(t)=\varphi_2(t,q(t),I(t)) \end{cases} \qquad \text{avec } \begin{cases} \varphi_1(t,q(t),I(t))=I(t) \\ \varphi_2(t,q(t),I(t))=-\frac{1}{LC}q(t) \end{cases} \)

Transformons le système en une EDO d’ordre 2: \( \begin{cases} L I'(t)+ \frac{1}{C}q(t)=E(t)\\ q'(t)=I(t)\implies q''(t)=I'(t) \implies Lq''(t)=LI'(t) \end{cases} \) donc \( L q''(t) =L I'(t) =-\frac{1}{C}q(t)+E(t) =-\frac{1}{C}q(t)+E(t) \) ainsi \( \boxed{Lq''(t)+\frac{1}{C}q(t)=0} \text{ et } \boxed{I(t)=q'(t)} \) avec \( \begin{cases} q(0)=q_0,\\ q'(0)=i_0. \end{cases} \)

On peut calculer la solution exacte à la main: \( \begin{cases} q''(t)+\frac{1}{LC}q(t)=0 \\ q(0)=q_0,\\ q'(0)=i_0 \end{cases} \) Il s’agit d’une EDO d’ordre 2 linéaire à coefficients constants et homogène. Le polynôme caractéristique est \(\lambda^2+\frac{1}{LC}=0\) et a les deux racines complexes conjugués \(\lambda_1=-i\frac{1}{\sqrt{LC}}\) et \(\lambda_2=i\frac{1}{\sqrt{LC}}\) ainsi la solution générale est \( q(t)= C_1 \cos\left( \frac{1}{\sqrt{LC}} t\right) + C_2 \sin\left(\frac{1}{\sqrt{LC}} t\right), \qquad C_1,C_2\in\mathbb{R}. \)

Sinon, on peut utiliser le module de calcul symbolique.

Code 
sp.var('t,t_0,q_0,i_0')
L = sp.Symbol('L',positive=True)
C = sp.Symbol('C',positive=True)

q = sp.Function('q')
edo = sp.Eq(sp.diff(q(t),t,t) , -q(t)/(L*C) )
display(edo)

solgen = sp.dsolve(edo,q(t)).simplify()
display(solgen)

\(\displaystyle \frac{d^{2}}{d t^{2}} q{\left(t \right)} = - \frac{q{\left(t \right)}}{C L}\)

\(\displaystyle q{\left(t \right)} = C_{1} \sin{\left(\frac{t}{\sqrt{C} \sqrt{L}} \right)} + C_{2} \cos{\left(\frac{t}{\sqrt{C} \sqrt{L}} \right)}\)

Avec les valeurs données:

Code 
pb = { L:2 , C:sp.S(1)/2 }
solgenLC = solgen.subs(pb)
display(solgenLC)

\(\displaystyle q{\left(t \right)} = C_{1} \sin{\left(t \right)} + C_{2} \cos{\left(t \right)}\)

Code 
t_0 = 0 
q_0 = 1
i_0 = 2

consts = sp.solve( [solgenLC.rhs.subs(t,t_0)-q_0 , solgenLC.rhs.diff(t).subs(t,t_0)- i_0],  dict=True)[0]
display(consts)

solparLC = solgenLC.subs(consts)
display(solparLC)

I = sp.Function('I')
I = solparLC.rhs.diff(t).simplify()
display(Math(r"I(t)="+sp.latex( I.subs(pb) )))
{C1: 2, C2: 1}

\(\displaystyle q{\left(t \right)} = 2 \sin{\left(t \right)} + \cos{\left(t \right)}\)

\(\displaystyle I(t)=- \sin{\left(t \right)} + 2 \cos{\left(t \right)}\)

On transforme les deux fonctions sympy en des fonctions numpy pour l’évaluation et l’affichage :

Code 
Qfunc = sp.lambdify(t,solparLC.rhs,'numpy')
Ifunc = sp.lambdify(t,I,'numpy')

tt = np.linspace(0,7,101)
yyQ = Qfunc(tt)
yyI = Ifunc(tt)

plt.figure(1,figsize=(17,12))

plt.subplot(2,3,1)
plt.plot(tt,yyQ);
plt.xlabel('$t$ (s)')
plt.ylabel('$q$ (C)')
plt.grid(True)

plt.subplot(2,3,2)
plt.plot(tt,yyI);
plt.xlabel('$t$ (s)')
plt.ylabel('$I$ (A)')
plt.grid(True)

plt.subplot(2,3,3)
plt.plot(yyQ,yyI);
plt.xlabel('$q$ (C)')
plt.ylabel('$I$ (A)')
plt.grid(True)

L = 2
C = 1/2

Ee=[0.5*(yyQ[i])**2/C for i in range(len(yyQ))]
plt.subplot(2,3,4)
plt.plot(tt,Ee);
plt.xlabel('$t$ (s)')
plt.ylabel('$E_C$ (?)')
plt.title("Electrique")
plt.grid(True)

Em=[0.5*L*(yyI[i])**2 for i in range(len(yyQ))]
plt.subplot(2,3,5)
plt.plot(tt,Em);
plt.xlabel('$t$ (s)')
plt.ylabel('$E_I$ (?)')
plt.title("Magnétique")
plt.grid(True)


Et=[Em[i]+Ee[i] for i in range(len(yyQ))]
plt.subplot(2,3,6)
plt.plot(tt,Et);
plt.xlabel('$t$ (s)')
plt.ylabel('$E_T$ (J)')
plt.title("Somme")
plt.grid(True)

Correction Q2

Code 
t0, tfinal = 0, 7
N = 51
q_0,i_0=1,2

phi1 = lambda t,q,qp : qp
phi2 = lambda t,q,qp : -q/(L*C)

tt = np. linspace(t0,tfinal,N)
yy1 = [Qfunc(t) for t in tt]
yy2 = [Ifunc(t) for t in tt]

\( \begin{cases} q_{n+1}=q_{n}+h\varphi_1(t_n,q_n,I_n) \\ I_{n+1}=I_n+h\varphi_2(t_n,q_n,I_n) \end{cases} \)

Code 
def EE(phi1,phi2,tt,q_0,i_0):
    h = tt[1]-tt[0]
    uu1 = [q_0]
    uu2 = [i_0]
    for i in range(len(tt)-1):
        uu1.append( uu1[i] + h*phi1(tt[i],uu1[i],uu2[i]) )
        uu2.append( uu2[i] + h*phi2(tt[i],uu1[i],uu2[i]) )
    return [uu1,uu2]
Code 
uu1E,uu2E = EE(phi1,phi2,tt,q_0,i_0)

plt.figure(1,figsize=(15,5))

plt.subplot(1,3,1)
plt.plot(tt,yy1,'r-',label=('Exacte'))
plt.plot(tt,uu1E,'b-o',label=('EE'))
plt.xlabel('$t$ (s)')
plt.ylabel('$q$ (C)')
plt.legend() 
plt.grid(True)

plt.subplot(1,3,2)
plt.plot(tt,yy2,'r-',label=('Exacte'))
plt.plot(tt,uu2E,'b-o',label=('EE'))
plt.xlabel('$t$ (s)')
plt.ylabel('$I$ (A)')
plt.legend() 
plt.grid(True)

plt.subplot(1,3,3)
plt.plot(yy1,yy2,'r-',label=('Exacte'))
plt.plot(uu1E,uu2E,'b-o',label=('EE'))
plt.xlabel('$q$ (C)')
plt.ylabel('$I$ (A)')
plt.legend() 
plt.grid(True)

Correction Q3
\( \begin{cases} \tilde q = q_{n}+h\varphi_1(t_n,q_n,I_n), \\ \tilde I = I_{n}+h\varphi_2(t_n,q_n,I_n), \\ q_{n+1}=q_{n}+\frac{h}{2}\Big(\varphi_1\big(t_n,q_n,I_n\big)+\varphi_1\big(t_{n+1},\tilde q,\tilde I\big)\Big) \\ I_{n+1}=I_n+\frac{h}{2}\Big(\varphi_2\big(t_n,q_n,I_n\big)+\varphi_2\big(t_{n+1},\tilde q,\tilde I\big)\Big) \end{cases} \)

Code 
def Heun(phi1,phi2,tt,q_0,i_0):
    h = tt[1]-tt[0]
    uu1 = [q_0]
    uu2 = [i_0]
    for i in range(len(tt)-1):
        uu1pred = uu1[i] + h*phi1(tt[i],uu1[i],uu2[i])
        uu2pred = uu2[i] + h*phi2(tt[i],uu1[i],uu2[i])
        uu1.append( uu1[i] + h/2*( phi1(tt[i],uu1[i],uu2[i])+phi1(tt[i],uu1pred,uu2pred) ) )
        uu2.append( uu2[i] + h/2*( phi2(tt[i],uu1[i],uu2[i])+phi2(tt[i],uu1pred,uu2pred) ) )
    return [uu1,uu2]
Code 
uu1H,uu2H = Heun(phi1,phi2,tt,q_0,i_0)

plt.figure(1,figsize=(15,5))

plt.subplot(1,3,1)
plt.plot(tt,yy1,'r-',label=('Exacte'))
plt.plot(tt,uu1H,'b-o',label=('H'))
plt.xlabel('$t$ (s)')
plt.ylabel('$q$ (C)')
plt.legend() 
plt.grid(True)

plt.subplot(1,3,2)
plt.plot(tt,yy2,'r-',label=('Exacte'))
plt.plot(tt,uu2H,'b-o',label=('H'))
plt.xlabel('$t$ (s)')
plt.ylabel('$I$ (A)')
plt.legend() 
plt.grid(True)

plt.subplot(1,3,3)
plt.plot(yy1,yy2,'r-',label=('Exacte'))
plt.plot(uu1H,uu2H,'b-o',label=('H'))
plt.xlabel('$q$ (C)')
plt.ylabel('$I$ (A)')
plt.legend() 
plt.grid(True)

Correction Q4

Code 
# Euler-Cromer, see fdm-book-4print.pdf page 47
# The scheme goes under several names: forward-backward scheme, semi-implicit Euler method, 
# semi-explicit Euler, symplectic Euler, Newton-Störmer-Verlet, and Euler-Cromer.
def EC(phi1,phi2,tt,q_0,i_0):
    h = tt[1]-tt[0]
    uu1 = [q_0]
    uu2 = [i_0]
    for i in range(len(tt)-1):
        uu2.append( uu2[i] + h*phi2(tt[i],uu1[i],uu2[i]) )
        uu1.append( uu1[i] + h*phi1(tt[i+1],uu1[i],uu2[i+1]) ) # en theorie uu1[i+1] mais n'intervient pas dans phi1
    return [uu1,uu2]
Code 
uu1EC,uu2EC = Heun(phi1,phi2,tt,q_0,i_0)

plt.figure(1,figsize=(15,5))

plt.subplot(1,3,1)
plt.plot(tt,yy1,'r-',label=('Exacte'))
plt.plot(tt,uu1EC,'b-o',label=('EC'))
plt.xlabel('$t$ (s)')
plt.ylabel('$q$ (C)')
plt.legend() 
plt.grid(True)

plt.subplot(1,3,2)
plt.plot(tt,yy2,'r-',label=('Exacte'))
plt.plot(tt,uu2EC,'b-o',label=('EC'))
plt.xlabel('$t$ (s)')
plt.ylabel('$I$ (A)')
plt.legend() 
plt.grid(True)

plt.subplot(1,3,3)
plt.plot(yy1,yy2,'r-',label=('Exacte'))
plt.plot(uu1EC,uu2EC,'b-o',label=('EC'))
plt.xlabel('$q$ (C)')
plt.ylabel('$I$ (A)')
plt.legend() 
plt.grid(True)