%reset -f
from IPython.display import display, Latex
from IPython.core.display import HTML
css_file = './custom.css'
HTML(open(css_file, "r").read())
Q1 [5 points]
Considérons le problème de Cauchy $$\begin{cases} y'(t)=y(t)-2e^{-t},\\ y(0)=a \end{cases}$$
- Calculer l'unique solution exacte à ce problème de Cauchy en fonction de $a$ pour $t>0$.
- Tracer la solution exacte sur l'intervalle $[0;4]$ pour $a=1$, $a=1.1$ et $a=0.9$.
- 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
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 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
%matplotlib inline
import sympy as sym
sym.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=0
u0=a
consts = sym.solve( [solgen.rhs.subs(t,t0)-u0 ], dict=True)[0]
print("Prise en compte de la condition initiale u({})={}:".format(t0,u0))
display(consts)
solpar=solgen.subs(consts).collect(t)
print("Solution du problème de Cauchy:")
display(solpar)
print("Affichage de trois solutions")
p = sym.plot(solpar.rhs.subs({a:1}), (t,0,4), line_color='b', legend = True, show=False)
p.extend(sym.plot(solpar.rhs.subs({a:1.1}), (t,0,4), line_color='m', legend = True, show=False))
p.extend(sym.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();
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.
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
Posons
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
import sympy as sym
sym.init_printing()
from IPython.display import display, Math
sym.var('t,t_0,q_0,i_0')
L=sym.Symbol('L',positive=True)
C=sym.Symbol('C',positive=True)
q=sym.Function('q')
edo=sym.Eq(L*sym.diff(q(t),t,t) , -1/C*q(t))
display(edo)
solgen=sym.dsolve(edo,q(t)).simplify()
display(solgen)
consts = sym.solve( [solgen.rhs.subs(t,t0)-q_0 , solgen.rhs.diff(t).subs(t,t_0)- i_0], dict=True)[0]
display(consts)
solpar=solgen.subs(consts)
display(solpar)
I= sym.Function('I')
I= solpar.rhs.diff(t).simplify()
display(Math(r"I(t)="+sym.latex( I )))
Avec les valeurs données:
pb={L:2,C:sym.S(1)/2}
solparLC=solpar.subs(pb)
display(solparLC)
init={t_0:0, q_0:1 , i_0:2}
solparLCinit=solparLC.subs(init)
display(solparLCinit)
display(Math(r"I(t)="+sym.latex( I.subs(init).subs(pb) )))
On transforme les deux fonctions sympy
en fonction vectorielle pour l'évaluation et l'affichage:
from matplotlib.pylab import *
Qfunc = sym.lambdify(t,solpar.rhs.subs(init).subs(pb),'numpy')
Ifunc = sym.lambdify(t,I.subs(init).subs(pb),'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=2
C=1/2
Ee=[0.5*(yyQ[i])**2/C for i in range(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])**2 for i in range(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 in range(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.
from matplotlib.pylab import *
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 = linspace(t0,tfinal,N)
yy1 = [Qfunc(t) for t in tt]
yy2 = [Ifunc(t) for t in tt]
Correction Q3
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]
uu1E,uu2E = EE(phi1,phi2,tt,q_0,i_0)
figure(1,figsize=(15,5))
subplot(1,3,1)
plot(tt,yy1,'r-',label=('Exacte'))
plot(tt,uu1E,'b-o',label=('EE'))
xlabel('$t$ (s)')
ylabel('$q$ (C)')
legend()
grid(True)
subplot(1,3,2)
plot(tt,yy2,'r-',label=('Exacte'))
plot(tt,uu2E,'b-o',label=('EE'))
xlabel('$t$ (s)')
ylabel('$I$ (A)')
legend()
grid(True)
subplot(1,3,3)
plot(yy1,yy2,'r-',label=('Exacte'))
plot(uu1E,uu2E,'b-o',label=('EE'))
xlabel('$q$ (C)')
ylabel('$I$ (A)')
legend()
grid(True)
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.
Correction Q4
$$
\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}
$$
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]
uu1H,uu2H = Heun(phi1,phi2,tt,q_0,i_0)
figure(1,figsize=(15,5))
subplot(1,3,1)
plot(tt,yy1,'r-',label=('Exacte'))
plot(tt,uu1H,'b-o',label=('H'))
xlabel('$t$ (s)')
ylabel('$q$ (C)')
legend()
grid(True)
subplot(1,3,2)
plot(tt,yy2,'r-',label=('Exacte'))
plot(tt,uu2H,'b-o',label=('H'))
xlabel('$t$ (s)')
ylabel('$I$ (A)')
legend()
grid(True)
subplot(1,3,3)
plot(yy1,yy2,'r-',label=('Exacte'))
plot(uu1H,uu2H,'b-o',label=('H'))
xlabel('$q$ (C)')
ylabel('$I$ (A)')
legend()
grid(True)
$$ \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} $$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, appellée Euler-Cromer ou Euler symplectique ou Newton-Störmer-Verlet, s'écrit comme suit:
Correction Q5
# 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]
uu1EC,uu2EC = Heun(phi1,phi2,tt,q_0,i_0)
figure(1,figsize=(15,5))
subplot(1,3,1)
plot(tt,yy1,'r-',label=('Exacte'))
plot(tt,uu1EC,'b-o',label=('EC'))
xlabel('$t$ (s)')
ylabel('$q$ (C)')
legend()
grid(True)
subplot(1,3,2)
plot(tt,yy2,'r-',label=('Exacte'))
plot(tt,uu2EC,'b-o',label=('EC'))
xlabel('$t$ (s)')
ylabel('$I$ (A)')
legend()
grid(True)
subplot(1,3,3)
plot(yy1,yy2,'r-',label=('Exacte'))
plot(uu1EC,uu2EC,'b-o',label=('EC'))
xlabel('$q$ (C)')
ylabel('$I$ (A)')
legend()
grid(True)