#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
TD5 au Lycée Lakanal (sujet rédigé par Arnaud Basson).
Traite de `la méthode d'Euler <https://fr.wikipedia.org/wiki/M%C3%A9thode_d%27Euler>`_ (:py:func:`tension`, :py:func:`tension_ordre2`) et du `pivot de Gauss <https://fr.wikipedia.org/wiki/%C3%89limination_de_Gauss-Jordan>`_ (cas tridiagonale, :py:func:`pivotTridiag`).
.. image:: exos_num__euler_ordre1_3.png
:align: center
:scale: 100%
- *Date :* vendredi 18 décembre 2015,
- *Auteur :* Lilian Besson, pour le cours d'informatique pour tous en prépa MP (http://perso.crans.org/besson/infoMP/),
- *Licence :* MIT Licence (http://lbesson.mit-license.org).
"""
from __future__ import print_function # Python 2 compatibility
import numpy as np
import matplotlib.pyplot as plt
# %% Exercice 1 : Circuits RC et RLC
if __name__ == '__main__':
print("\n\nExercice 1 : Circuits RC et RLC")
# Question 1
[docs]def tension(h, tmax, ue, R, C):
r"""
Résoud numériquement l'équation différentielle :eq:`circuit_rc` par `la méthode d'Euler <https://fr.wikipedia.org/wiki/M%C3%A9thode_d%27Euler>`_ pour :math:`t \in [0, t_{max}]` avec un pas de temps h.
.. math:: R C \frac{d u_c}{d t}(t) + u_c(t) = u_e(t)
:label: circuit_rc
Pour utiliser un schéma d'Euler, il nous faut écrire cette équation sous la forme :math:`u'(t) = f(t, u)` :
.. math::
\frac{d u_c}{d t}(t) &= (u_e(t) - u_c(t)) / (R \times C) \\
f(t, u_c(t)) &= (u_e(t) - u_c(t)) / (R \times C)
Et dès lors, le schéma de mise à jour d'Euler s'écrit :
.. math::
u_{i+1} = u_i + (t_{i+1} - t_i) \times f(t_i, u_i)
- Renvoit un tableau (``list``) contenant les valeurs de :math:`u_c` calculées aux instants :math:`t_i = i \times h` (pour :math:`0 \leq i \leq t_{max} / h)`.
- *Arguments* :
- ``h`` est un un pas de temps (:math:`h > 0`),
- ``tmax`` est la durée totale de la simulation numérique (:math:`t_{max} > 0`),
- et ``ue`` est *une fonction* :math:`u_e t \mapsto u_e(t)`, est la tension d'entrée (tension aux bornes de l'ensemble résistance + condensateur).
- *Hypothèse* : On résoudra l'équation pour :math:`t > 0` en supposant qu'à :math:`t = 0` le condensateur est déchargé (:math:`u_c(0) = 0`).
- Exemple (question 2.a), avec une tension d'entrée *sinusoïdale* :
>>> R = 1e6; C = 10e-9
>>> tmax = 0.2
>>> h = tmax / 1000 # On choisit 1000 points
>>> uemax = 10; freq = 50
>>> ue = lambda t: uemax * np.sin(t * freq)
>>> ue(tmax)
-5.4402111088936973
>>> uc = tension(h, tmax, ue, R, C)
>>> uc[-1]
-0.89989324797564696
.. image:: exos_num__euler_ordre1_1.png
:align: center
:scale: 100%
- Exemple (question 2.b), avec une tension d'entrée *triangulaire*.
.. image:: exos_num__euler_ordre1_2.png
:align: center
:scale: 100%
- Exemple (question 2.c), avec une tension d'entrée *créneaux* :
>>> R = 1e6; C = 10e-9
>>> tmax = 0.2
>>> h = tmax / 1000 # On choisit 1000 points
>>> uemax = 10; freq = 50
>>> ue = lambda t: uemax * np.sin(1e-12 + t * freq) / np.abs(np.sin(1e-12 + t * freq))
>>> ue(tmax)
-10.0
>>> uc = tension(h, tmax, ue, R, C)
>>> uc[-1]
-3.5593830469494541
.. image:: exos_num__euler_ordre1_3.png
:align: center
:scale: 100%
.. note::
Pour plus de détails sur les circuits RC, voir `cette page <https://fr.wikipedia.org/wiki/Oscillateur_harmonique#Oscillateurs_.C3.A9lectriques>`_.
"""
i_max = int(tmax / h)
# On créé un tableau qu'on va remplir petit à petit
uc = [0 for i in range(0, i_max)]
# f = lambda i: (ue[i] - uc[i]) / (R * C)
for i in range(0, i_max-1):
# On applique la relation d'Euler
uc[i+1] = uc[i] + h * (ue(h*i) - uc[i]) / (R * C)
# y_i+1 = y_i + (x_i+1 - x_i) * f(x_i, y_i)
# f(x_i, y_i) = f(i) = (ue(t[i]) - uc[i]) / (R * C)
return uc
# Question 2 : Simulations
if False and __name__ == '__main__':
R = 1e6
C = 10e-9
tmax = 0.2
h = tmax / 1000 # On choisit 1000 points
uemax = 10
freq = 50
# 2.a) ue est une sinusoïdale
def ue(t):
""" Tension sinusoïdale."""
return uemax * np.sin(t * freq)
uc = tension(h, tmax, ue, R, C)
t = np.arange(0, tmax, h)
plt.figure()
plt.plot(t, [ue(ti) for ti in t], label="Tension d'entrée sinusoïdale $u_e(t)$")
plt.plot(t, uc, label="Tension du condensateur $u_c(t)$")
plt.legend()
plt.xlabel(r"Temps $0 \leq t \leq {}$".format(tmax))
plt.xlim([0, tmax*1.01])
plt.ylabel(r"Tension du condensateur $u_c(t)$")
plt.ylim([-uemax*1.05, uemax*1.05])
plt.title(r"2.a - Résolution d'une équation différentielle d'ordre 1 via la méthode d'Euler")
plt.show()
# 2.b) ue est une triangulaire
def ue(t):
""" Tension triangulaire."""
return uemax * (-1 + 4*abs(((t*freq/(2*np.pi)))%1 - 0.5))
uc = tension(h, tmax, ue, R, C)
t = np.arange(0, tmax, h)
plt.figure()
plt.plot(t, [ue(ti) for ti in t], label="Tension d'entrée triangulaire $u_e(t)$")
plt.plot(t, uc, label="Tension du condensateur $u_c(t)$")
plt.legend()
plt.xlabel(r"Temps $0 \leq t \leq {}$".format(tmax))
plt.xlim([0, tmax*1.01])
plt.ylabel(r"Tension du condensateur $u_c(t)$")
plt.ylim([-uemax*1.05, uemax*1.05])
plt.title(r"2.b - Résolution d'une équation différentielle d'ordre 1 via la méthode d'Euler")
plt.show()
# 2.c) ue est une créneaux
def ue(t):
""" Tension créneaux."""
return uemax * (-1) ** (4*(t*freq/(2*np.pi)) // 2)
uc = tension(h, tmax, ue, R, C)
t = np.arange(0, tmax, h)
plt.figure()
plt.plot(t, [ue(ti) for ti in t], label="Tension d'entrée créneaux $u_e(t)$")
plt.plot(t, uc, label="Tension du condensateur $u_c(t)$")
plt.legend()
plt.xlabel(r"Temps $0 \leq t \leq {}$".format(tmax))
plt.xlim([0, tmax*1.01])
plt.ylabel(r"Tension du condensateur $u_c(t)$")
plt.ylim([-uemax*1.05, uemax*1.05])
plt.title(r"2.c - Résolution d'une équation différentielle d'ordre 1 via la méthode d'Euler")
plt.show()
# Question 3
# Même genre de fonction, même genre de graphiques...
[docs]def tension_ordre2(h, tmax, ue, R, L, C):
r"""
Résoud numériquement l'équation différentielle d'ordre 2 :eq:`circuit_lrc` par `la méthode d'Euler <https://fr.wikipedia.org/wiki/M%C3%A9thode_d%27Euler>`_ pour :math:`t \in [0, t_{max}]` avec un pas de temps h.
.. math::
:label: circuit_lrc
L C \frac{d^2 u_c}{d t^2}(t) + R C \frac{d u_c}{d t}(t) + u_c(t) = u_e(t)
Pour utiliser un schéma d'Euler, il nous faut écrire cette équation (linéaire d'ordre 2) sous la forme :math:`U'(t) = f(t, U)`, où :math:`U = [u, u']` :
.. math::
U(t) &= [ u_c(t), \frac{d u_c}{d t}(t) ] \\
U(t) &= \Bigl[\begin{smallmatrix} u_c(t) \\ \frac{d u_c}{d t}(t) \end{smallmatrix} \Bigr] \;\;\text{(comme un vecteur)} \\
\frac{d U}{d t}(t) &= \Bigl[ U[2](t), (u_e(t) - U[1](t))/(L\times C) - R/L \times U[2](t) \Bigr] \\
F(t, U) &= \Bigl[ U[2](t), (u_e(t) - U[1](t))/(L \times C) - R/L \times U[2](t) \Bigr]
Et dès lors, le schéma de mise à jour d'Euler s'écrit :
.. math::
U_{i+1} = U_i + (t_{i+1} - t_i) \times F(t_i, U_i)
- Cette fois l'opérateur :math:`F: (t, U) \mapsto f(t, U)` est vectoriel (de taille :math:`2 \times 2`), et affine de qui plus est (mais ça n'a aucune importance !).
- Renvoit un tableau (``list``) contenant les valeurs de :math:`u_c` calculées aux instants :math:`t_i = i \times h` (pour :math:`0 \leq i \leq t_{max} / h)`.
- *Arguments* :
- ``h`` est un un pas de temps (:math:`h > 0`),
- ``tmax`` est la durée totale de la simulation numérique (:math:`t_{max} > 0`),
- et ``ue`` est *une fonction* :math:`u_e t \mapsto u_e(t)`, est la tension d'entrée (tension aux bornes de l'ensemble résistance + condensateur).
- *Hypothèse* : On résoudra l'équation pour :math:`t > 0` en supposant qu'à :math:`t = 0` le condensateur est déchargé (:math:`u_c(0) = 0`) et que :math:`\frac{d u_c}{d t}(t=0) = 0` aussi (sinon il nous faut une condition numérique précise).
- Exemple (question 3.a), avec une tension d'entrée *sinusoïdale* :
>>> R = 1e3; L = 1; C = 100e-6
>>> tmax = 0.2
>>> h = tmax / 1000 # On choisit 1000 points
>>> uemax = 10; freq = 50
>>> ue = lambda t: uemax * np.sin(t * freq)
>>> ue(tmax)
-5.4402111088936973
>>> # Pour un circuit RC (ordre 1) :
>>> uc = tension(h, tmax, ue, R, C)
>>> uc[-1]
1.686614450937342
>>> # Pour un circuit RLC (ordre 2) :
>>> uc = tension_ordre2(h, tmax, ue, R, L, C)
>>> uc[-1] # Moins d'atténuation
1.7632486787067487
.. image:: exos_num__euler_ordre2_1.png
:align: center
:scale: 100%
- Exemple (question 3.b), avec une tension d'entrée *triangulaire*.
.. image:: exos_num__euler_ordre2_2.png
:align: center
:scale: 100%
- Exemple (question 3.c), avec une tension d'entrée *créneaux* :
>>> R = 1e3; L = 1; C = 100e-6
>>> tmax = 0.2
>>> h = tmax / 1000 # On choisit 1000 points
>>> uemax = 10; freq = 50
>>> ue = lambda t: uemax * np.sin(1e-12 + t * freq) / np.abs(np.sin(1e-12 + t * freq))
>>> ue(tmax)
-10.0
>>> # Pour un circuit RC (ordre 1) :
>>> uc = tension(h, tmax, ue, R, C)
>>> uc[-1]
2.0756606150715764
>>> # Pour un circuit RLC (ordre 2) :
>>> uc = tension_ordre2(h, tmax, ue, R, L, C)
>>> uc[-1] # Moins d'atténuation
2.2083917041910097
.. image:: exos_num__euler_ordre2_3.png
:align: center
:scale: 100%
.. note::
Pour plus de détails sur les circuits RLC, voir `cette page <https://fr.wikipedia.org/wiki/Oscillateur_harmonique#Oscillateurs_.C3.A9lectriques>`_.
"""
i_max = int(tmax / h)
# On créé un tableau qu'on va remplir petit à petit
uc = np.zeros((i_max, 2))
for i in range(0, i_max-1):
# On applique la relation d'Euler vectorielle
# [U_2(t), (ue - U_1(t))/LC - R/L U_2(t)]
uc[i+1, 0] = uc[i, 0] + h * uc[i, 1]
uc[i+1, 1] = uc[i, 1] + h * (((ue(i*h) - uc[i, 0]) / (L * C)) - ((R/L) * uc[i, 1]))
# On ne renvoit que les positions pas les vitesses
return uc[:, 0]
# Question 3 : Simulations
if False and __name__ == '__main__':
R = 1e3
L = 1
C = 100e-6
tmax = 0.2
h = tmax / 1000 # On choisit 1000 points
uemax = 10
freq = 50
# 2.a) ue est une sinusoïdale
def ue(t):
""" Tension sinusoïdale."""
return uemax * np.sin(t * freq)
uc = tension_ordre2(h, tmax, ue, R, L, C)
t = np.arange(0, tmax, h)
plt.figure()
plt.plot(t, [ue(ti) for ti in t], label="Tension d'entrée sinusoïdale $u_e(t)$")
plt.plot(t, uc, label="Tension du condensateur $u_c(t)$")
plt.legend()
plt.xlabel(r"Temps $0 \leq t \leq {}$".format(tmax))
plt.xlim([0, tmax*1.01])
plt.ylabel(r"Tension du condensateur $u_c(t)$")
plt.ylim([-uemax*1.05, uemax*1.05])
plt.title(r"3.a - Résolution d'une équation différentielle d'ordre 2 via la méthode d'Euler")
plt.show()
# 2.b) ue est une triangulaire
def ue(t):
""" Tension triangulaire."""
return uemax * (-1 + 4*abs(((t*freq/(2*np.pi)))%1 - 0.5))
uc = tension_ordre2(h, tmax, ue, R, L, C)
t = np.arange(0, tmax, h)
plt.figure()
plt.plot(t, [ue(ti) for ti in t], label="Tension d'entrée triangulaire $u_e(t)$")
plt.plot(t, uc, label="Tension du condensateur $u_c(t)$")
plt.legend()
plt.xlabel(r"Temps $0 \leq t \leq {}$".format(tmax))
plt.xlim([0, tmax*1.01])
plt.ylabel(r"Tension du condensateur $u_c(t)$")
plt.ylim([-uemax*1.05, uemax*1.05])
plt.title(r"3.b - Résolution d'une équation différentielle d'ordre 2 via la méthode d'Euler")
plt.show()
# 2.c) ue est une créneaux
def ue(t):
""" Tension créneaux."""
return uemax * (-1) ** (4*(t*freq/(2*np.pi)) // 2)
uc = tension_ordre2(h, tmax, ue, R, L, C)
t = np.arange(0, tmax, h)
plt.figure()
plt.plot(t, [ue(ti) for ti in t], label="Tension d'entrée créneaux $u_e(t)$")
plt.plot(t, uc, label="Tension du condensateur $u_c(t)$")
plt.legend()
plt.xlabel(r"Temps $0 \leq t \leq {}$".format(tmax))
plt.xlim([0, tmax*1.01])
plt.ylabel(r"Tension du condensateur $u_c(t)$")
plt.ylim([-uemax*1.05, uemax*1.05])
plt.title(r"3.c - Résolution d'une équation différentielle d'ordre 2 via la méthode d'Euler")
plt.show()
# %% Exercice 2 : Le pivot de Gauss pour une matrice tridiagonale
if __name__ == '__main__':
print("\n\nExercice 2 : Le pivot de Gauss pour une matrice tridiagonale")
from numpy.linalg import inv, solve
# Question a
[docs]def pivotTridiag(A, b, verb=False):
r"""
Calcule la solution du système linéaire :math:`Ax = b`, en appliquant l'algorithme `du pivot de Gauss <https://fr.wikipedia.org/wiki/%C3%89limination_de_Gauss-Jordan>`_, qui se simplifie largement dans le cas tridiagonal.
- *Hypothèses* : on suppose que A est bien inversible, et qu'aucun échange de ligne n'est nécessaire au cours de la résolution du système linéaire.
- *Arguments* :
- ``A`` : matrice tridiagonale, de taille :math:`n \times n`, supposée inversible,
- ``b`` : un vecteur b de taille :math:`n` (ie. :math:`b \in \mathbb{R}^n`).
- Complexité *temporelle* : :math:`O(n)`,
- Complexité *mémoire* : :math:`O(n^2)`.
- Exemple :
>>> A = np.array([[ 1, -2, 0], [ 3, 4, -5], [ 0, 6, 7]])
>>> b = np.array([-1, -1, 0])
>>> x = pivotTridiag(A, b)
>>> x
array([-0.72, 0.14, -0.12])
"""
# A est une matrice tridiagonale inversible et b un vecteur
n = len(A)
# 1e étape: mise sous forme triangulaire
for k in range(n-1):
pivot = A[k, k]
A[k, k] = 1.0 # normalisation
A[k, k+1] /= pivot
b[k] /= pivot
c = A[k+1, k]
# L[k+1] <-- L[k+1] - a[k+1, k]*L[k]
A[k+1, k+1] -= c*A[k, k+1]
A[k+1, k] = 0.0
b[k+1] -= c*b[k]
# Normalisation du dernier pivot
b[n-1] /= A[n-1, n-1]
A[n-1, n-1] = 1.0
# 2e étape: calcul de x (le système est triangulaire supérieur avec des 1 sur la diagonale)
x = np.zeros(n)
x[n-1] = b[n-1]
for k in range(n-2, -1, -1):
x[k] = b[k] - A[k, k+1]*x[k+1]
return x
# Question b
[docs]def lireMatriceA(nomFichier):
r"""
Lit la matrice A stockée dans le fichier ``nomFichier``, en suivant le format décrit dans l'énoncé.
- Les coefficients de la matrice :math:`A` sont stockés dans un fichier texte ``nomFichier`` (e.g. ``matrice.txt``).
- La première ligne du fichier contient la valeur de :math:`n` et les :math:`3n - 2` lignes suivantes contiennent les valeurs des :math:`a_k` , puis des :math:`b_k` puis des :math:`c_k` , à raison d’un seul coefficient par ligne.
Interlude sur la lecture et écriture de fichiers texte en python :
- On créé deux fichiers inutiles pour tester leur lecture et écriture (*juste une démo*).
.. runblock:: console
$ echo "Fichier en lecture" > /tmp/read.txt
$ echo "Fichier en ecriture (ecrase)" > /tmp/write.txt
- On rappelle qu'il faut utiliser la fonction :py:func:`open` pour ouvrir un fichier, en mode ``"w"`` pour écrire (*write*) ou ``"r"`` pour lire (*read*) :
>>> fichierEcriture = open("/tmp/write.txt", "w")
>>> fichierEcriture.write("Ecrase!") # Renvoit le nombre de caracteres ecrits
7
>>> fichierLecture = open("/tmp/read.txt", "r")
>>> print(fichierLecture.readline())
OK
<BLANKLINE>
.. runblock:: console
$ echo cat /tmp/write.txt # Le contenu a change
Retour sur l'écriture d'une matrice :
- Exemple : le fichier `exos_num_matrice_A.txt <../exos_num_matrice_A.txt>`_ contient un matrice tridiagonale :math:`A` de taille :math:`3 \times 3` :
.. math:: A = \Bigl[\begin{smallmatrix} 1 & -2 & 0 \\ 3 & 4 & -5 \\ 0 & 6 & 7 \end{smallmatrix} \Bigr]
.. runblock:: console
$ echo "Contenu du fichier 'exos_num_matrice_A.txt' :"
$ cat exos_num_matrice_A.txt
- On peut donc lire le contenu du fichier et le sauvegarder dans la matrice (numpy array) A :
>>> A = lireMatriceA("exos_num_matrice_A.txt")
>>> A
array([[ 1., -2., 0.],
[ 3., 4., -5.],
[ 0., 6., 7.]])
.. note::
Plus généralement, numpy propose les fonctions :py:func:`numpy.loadtxt` (`numpy.loadtxt (doc) <https://docs.scipy.org/doc/numpy/reference/generated/numpy.loadtxt.html>`_) et :py:func:`numpy.savetxt` (`numpy.savetxt (doc) <https://docs.scipy.org/doc/numpy/reference/generated/numpy.savetxt.html>`_) pour lire et sauvegarder des ``array`` dans un fichier texte (extension ``.txt``).
Et un format `encore plus compressé <https://docs.scipy.org/doc/numpy/neps/npy-format.html>`_ (espace réduit, sauvegarde et lecture plus rapides) est un format binaire spécialisé (extension ``.npy``), utilisable avec les fonctions :py:func:`numpy.load` (`numpy.load (doc) <https://docs.scipy.org/doc/numpy/reference/generated/numpy.load.html>`_) et :py:func:`numpy.save` (`numpy.save (doc) <https://docs.scipy.org/doc/numpy/reference/generated/numpy.save.html>`_).
Voir `cette page pour plus de détails <https://docs.scipy.org/doc/numpy/neps/npy-format.html#use-cases>`_ ou `cette page pour d'autres méthodes <https://scipy.github.io/old-wiki/pages/Cookbook/InputOutput#numPy>`_ (en anglais).
*Note :* ces notions plus avancées sont hors programme, bien sûr.
"""
# Ouvrir en mode lecture ("r")
f = open(nomFichier, "r")
n = int(f.readline()[:-1])
A = np.zeros((n, n))
A[0, 0] = float(f.readline()[:-1])
for i in range(1, n):
A[i-1, i] = float(f.readline()[:-1]) # b_i
A[i, i-1] = float(f.readline()[:-1]) # c_i
A[i, i] = float(f.readline()[:-1]) # a_i+1
f.close()
return A
[docs]def enregistrerSolution(A, b, nomFichier):
"""
Prenant en arguments une matrice tridiagonale inversible A, un vecteur b et le nom d'un fichier, qui calcule la solution du système Ax = b puis écrit dans le fichier indiqué les valeurs des coefficients de x (un coefficient par ligne).
- *Attention* : le contenu du fichier ``nomFichier`` sera écrasé !
- Exemple : résolution et sauvegarde de la solution d'un système linéaire de taille :math:`3` :
>>> import numpy as np
>>> A = lireMatriceA("exos_num_matrice_A.txt")
>>> A
array([[ 1., -2., 0.],
[ 3., 4., -5.],
[ 0., 6., 7.]])
>>> b = np.array([-1, -1, 0])
>>> x = pivotTridiag(A, b)
>>> x
array([-0.72, 0.14, -0.12])
>>> enregistrerSolution(A, b, 'exos_num_vecteur_x.txt')
Le vecteur x = [-0.72 0.14 -0.12] a bien ete sauvegarde dans le fichier 'exos_num_vecteur_x.txt'.
- Le fichier `exos_num_vecteur_x.txt <../exos_num_vecteur_x.txt>`_ contient désormais un vecteur de taille :math:`3` :
.. runblock:: console
$ echo "Contenu du fichier 'exos_num_vecteur_x.txt' :"
$ cat exos_num_vecteur_x.txt
"""
x = pivotTridiag(A, b)
n = np.size(x) # Ou np.shape(x)[0]
# Ouvrir en mode écriture ("w")
f = open(nomFichier, "w")
f.write(str(n) + '\n')
for i in range(n):
f.write(str(x[i]) + '\n')
f.close()
print("Le vecteur x = {} a bien ete sauvegarde dans le fichier '{}'.".format(x, nomFichier))
# %% Fin de la correction.
if __name__ == '__main__':
from doctest import testmod
print("Test automatique de toutes les doctests ecrites dans la documentation (docstring) de chaque fonction :")
testmod(verbose=True)
testmod()
# testmod()
print("\nPlus de details sur ces doctests peut etre trouve dans la documentation de Python:\nhttps://docs.python.org/3/library/doctest.html (en anglais)")
# Fin de exos_num.py