Q1 [5 points]
Considérons le problème de Cauchy \(\begin{cases}
y'(t)=y(t)-2e^{-t},\\
y(0)=a
\end{cases}\) 1. Calculer l’unique solution exacte à ce problème de Cauchy en fonction de \(a\) pour \(t>0\). 2. Tracer la solution exacte 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 Q1
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}\),
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 mais attention, pour ne pas avoir des mauvaises intéractions entre le module matplotlib et le module sympy, on utilisera l’import de sympy avec un alias: import sympy as sym.
Code
%reset -f%matplotlib inlineimport sympy as symsym.init_printing()sym.var('t,a')u = sym.Function('u')f = u(t)-2*sym.exp(-t)edo = sym.Eq(sym.diff(u(t),t),f)print("EDO:")display(edo)solgen = sym.dsolve(edo,u(t)).simplify()print("Solution générale:")display(solgen)t0 =0u0 = a consts = sym.solve( [solgen.rhs.subs(t,t0)-u0 ], dict=True)[0]print(f"Prise en compte de la condition initiale u({</span>t0<span class="sc">})={</span>u0<span class="sc">}:")display(consts)solpar = solgen.subs(consts)print("Solution du problème de Cauchy:")display(solpar)print("Affichage de trois solutions")p = sym.plot(solpar.rhs.subs({a:<span class="dv">1</span>}), (t,0,4), line_color='b', legend =True, show=False)p.extend(sym.plot(solpar.rhs.subs({a:<span class="fl">1.1</span>}), (t,0,4), line_color='m', legend =True, show=False))p.extend(sym.plot(solpar.rhs.subs({a:<span class="fl">0.9</span>}), (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:
Solution générale:
Prise en compte de la condition initiale u(0)=a:
Solution du problème de Cauchy:
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.
1.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.
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).
Q2 [3 points]
É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 ensuite la solution exacte.
Afficher \(t\mapsto q(t)\), \(t\mapsto i(t)\) et \(q\mapsto i(q)\).
Correction Q2
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 mais attention, pour ne pas avoir des mauvaises intéractions entre le module matplotlib et le module sympy, on utilisera l’import de sympy avec un alias: import sympy as sym.
On transforme les deux fonctions sympy en fonction vectorielle pour l’évaluation et l’affichage:
Code
from matplotlib.pylab import*Qfunc = sym.lambdify(t,solparLC.rhs,'numpy')Ifunc = sym.lambdify(t,I,'numpy')tt = linspace(0,7,101)yyQ = Qfunc(tt)yyI = Ifunc(tt)figure(1,figsize=(17,12))subplot(2,3,1)plot(tt,yyQ);xlabel('$t$ (s)')ylabel('$q$ (C)')grid(True)subplot(2,3,2)plot(tt,yyI);xlabel('$t$ (s)')ylabel('$I$ (A)')grid(True)subplot(2,3,3)plot(yyQ,yyI);xlabel('$q$ (C)')ylabel('$I$ (A)')grid(True)L =2C =1/2Ee=[0.5*(yyQ[i])**2/C for i inrange(len(yyQ))]subplot(2,3,4)plot(tt,Ee);xlabel('$t$ (s)')ylabel('$E_C$ (?)')title("Electrique")grid(True)Em=[0.5*L*(yyI[i])**2for i inrange(len(yyQ))]subplot(2,3,5)plot(tt,Em);xlabel('$t$ (s)')ylabel('$E_I$ (?)')title("Magnétique")grid(True)Et=[Em[i]+Ee[i] for i inrange(len(yyQ))]subplot(2,3,6)plot(tt,Et);xlabel('$q$ (C)')ylabel('$E_T$ (A)')title("Somme")grid(True)
Q3 [4 points]
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.
Code
from matplotlib.pylab import*t0, tfinal =0, 7N =51q_0,i_0=1,2phi1 =lambda t,q,qp : qpphi2 =lambda t,q,qp : -q/(L*C)tt = linspace(t0,tfinal,N)yy1 = [Qfunc(t) for t in tt]yy2 = [Ifunc(t) for t in tt]
Q4 [4 points]
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 de Heun pour \(t\in[0;7]\) avec \(N=101\) points.
Q5 [4 points] 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 de Euler-Cromer pour \(t\in[0;7]\) pour \(N=101\) points.
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 Q5
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 inrange(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 phi1return [uu1,uu2]
