Code source de exos_num

#!/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