In [2]:
%reset -f
from IPython.display import display, Latex
from IPython.core.display import HTML
css_file = './custom.css'
HTML(open(css_file, "r").read())
Out[2]:

M62 : Correction Examen 2019 Session 2

Choix d'un schéma

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}$,

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
In [19]:
%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();

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:
$\displaystyle \left\{ C_{1} : a - 1\right\}$ 
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.

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. Tuned_circuit_animation_3.gif

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

Posons

  • $L=2$H,
  • $C=0.5$F,
  • $q(0)=1$C,
  • $i(0)=2$A.

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
In [20]:
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 )))
$\displaystyle L \frac{d^{2}}{d t^{2}} q{\left(t \right)} = - \frac{q{\left(t \right)}}{C}$
$\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)}$
$\displaystyle \left\{ C_{1} : \frac{\sqrt{C} \sqrt{L} i_{0}}{\cos{\left(\frac{t_{0}}{\sqrt{C} \sqrt{L}} \right)}} + q_{0} \tan{\left(\frac{t_{0}}{\sqrt{C} \sqrt{L}} \right)}, \ C_{2} : q_{0}\right\}$
$\displaystyle q{\left(t \right)} = q_{0} \cos{\left(\frac{t}{\sqrt{C} \sqrt{L}} \right)} + \left(\frac{\sqrt{C} \sqrt{L} i_{0}}{\cos{\left(\frac{t_{0}}{\sqrt{C} \sqrt{L}} \right)}} + q_{0} \tan{\left(\frac{t_{0}}{\sqrt{C} \sqrt{L}} \right)}\right) \sin{\left(\frac{t}{\sqrt{C} \sqrt{L}} \right)}$
$\displaystyle I(t)=\frac{\sqrt{C} \sqrt{L} i_{0} \cos{\left(\frac{t}{\sqrt{C} \sqrt{L}} \right)} - q_{0} \sin{\left(\frac{t - t_{0}}{\sqrt{C} \sqrt{L}} \right)}}{\sqrt{C} \sqrt{L} \cos{\left(\frac{t_{0}}{\sqrt{C} \sqrt{L}} \right)}}$

Avec les valeurs données:

In [21]:
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) )))
$\displaystyle q{\left(t \right)} = q_{0} \cos{\left(t \right)} + \left(\frac{i_{0}}{\cos{\left(t_{0} \right)}} + q_{0} \tan{\left(t_{0} \right)}\right) \sin{\left(t \right)}$
$\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 fonction vectorielle pour l'évaluation et l'affichage:

In [22]:
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.

In [23]:
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

$$ \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} $$
In [24]:
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]
In [25]:
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} $$

In [26]:
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]
In [27]:
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)

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:

$$ \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

In [28]:
# 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]
In [29]:
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)