Équation de la chaleur

In [ ]:
%pylab inline
Populating the interactive namespace from numpy and matplotlib

Remarque préliminaire sur l'utilisation de symboles exotiques dans les noms de variables Python

Un programme Python est un texte ordinaire, c'est-à-dire une suite de caractères Unicode. Ce texte doit bien entendu se conformer à la syntaxe Python. Cette syntaxe spécifie en particulier les caractères autorisés pour les noms des variables. En fait, l'immense majorité des caractères sont autorisés. Les symboles interdits sont ceux qui ont une signification particulière pour Python (par exemple les symboles de ponctuation ou les symboles des opérations arithmétiques).

In [ ]:
उत्पादन = 5
print(उत्पादन * 2)
10

L'utilisation de symboles régionaux pose plusieurs problèmes :

  • L'ordinateur et le système d'exploitation utilisés par le lecteur du programme ne disposent pas des polices de caractères permettant de visualiser correctement le texte.
  • L'ordinateur et le système d'exploitation ne permettent pas de saisir facilement ces caractères régionaux.
  • La personne qui lit le programme n'est pas familière avec ces caractères.

Pour ces raisons, l'usage des professionnels est de se restreindre aux caractères ASCII pour les noms des variables. Ceci exclue donc aussi les lettres latines accentuées.

Lors des épreuves écrites des concours par contre, il n'y a a priori pas d'objection au fait d'utiliser des accents ou des lettres grecs sur la copie. Toutefois, il peut se trouver des correcteurs un peu sourcilleux qui pourraient vous le reprocher. Je vous conseille donc de n'utiliser que des caractères ASCII sur votre copie pendant les concours.

Puisque nous venons d'énoncer un principe, nous allons commencer par y faire une exception en utilisant des lettres grecs. Les problèmes mentionnés ci-dessus n'en sont pas vraiment dans ce TP puisque le code que nous allons écrire n'est destiné à n'être utilisé que dans ce notebook. L'intérêt est d'avoir des noms de variable python qui collent au plus près avec les notations physiques utilisées par ailleurs dans le texte du notebook.

Pour saisir une lettre grec $\alpha$ dans le notebook, tapez \alpha suivi d'une tabulation :

In [ ]:
α = 1

Le problème physique

On considère le problème bidimensionnel de la conduction de la chaleur dans une plaque carrée de 1 mètre de longueur.

Modèlisation du problème

Modéle continu

Le champ de température $u(t, x, y)$ vérifie l'équation de la chaleur :

$$ \frac{\partial u}{\partial t} = \alpha\left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right) \quad\quad (t, x, y) \in [0,+\infty[\times[0, 1]^2 $$

La température initiale de la plaque est donnée par la fonction $u_{\text{init}}(x, y)$ :

$$ u(0, x, y) = u_{\text{init}}(x, y)\quad\quad (x, y) \in \times[0, 1]^2 $$

La température sur les bords de la plaque est maintenue constante en fonction du temps. Elle est définie par les fonctions $u_{\text{left}}(y)$, $u_{\text{right}}(y)$, $u_{\text{bottom}}(x)$ et $u_{\text{top}}(x)$ :

$$ u(t, 0, y) = u_{\text{left}}(y)\quad \text{et}\quad u(t, 1, y) = u_{\text{right}}(y)\quad\quad (t, y) \in [0,+\infty[\times[0, 1] \\ u(t, x, 0) = u_{\text{bottom}}(x)\quad \text{et}\quad u(t, y, 1) = u_{\text{top}}(x)\quad\quad (t, x) \in [0,+\infty[\times[0, 1] $$

Modèle discret

Les ordinateurs ne connaissent que le fini et le discret. Pour résoudre informatiquement ce problème il est nécessaire de remplacer la modélisation continue ci-dessus par une modélisation discrète.

Une façon de discrétiser le problème est de représenter la plaque par un maillage carré de $n+1$ points sur $n+1$ points. On note $\Delta h$ le pas de discrétisation. On a donc $n\Delta h=1$. On note $x_i=i\Delta h$ pour $0\leq i\leq n$ et $y_j=j\Delta h$ pour $0\leq j\leq n$.

On discrétise également le temps en un nombre fini d'instants d'un intervalle borné $[0, t_{\text{max}}]$. On note $\Delta t$ le pas de temps et $t_k=k\Delta t$ pour $k=0,1,2,\dots$ les instants successifs de la discrétisation. On note $U_{i,j}^k$ la température de la plaque au point $(x_i,\,y_j)$ à l'instant $t_k$.

Discrétisation de l'équation aux dérivées partielles

Il s'agit maintenant de discrétiser l'EDP afin d'établir des équations vérifiées par les $U_{i,j}^k$. Par définition de la dérivée on a

$$\frac{\partial u}{\partial t}(t,x,y) = \lim_{\Delta t\to 0} \frac{u(t+\Delta t,x,y) - u(t,x,y)}{\Delta t}$$

On peut donc estimer que $$\frac{\partial u}{\partial t}(t_k,x_i,y_j)\approx\frac{U_{i,j}^{k+1}-U_{i,j}^k}{\Delta t}$$

Passons aux dérivées par rapport aux variables d'espace. La formule de Taylor à l'ordre 3 donne

$$ u(t,x+\Delta x,y) = u(t,x,y) + \Delta x\frac{\partial u}{\partial x}(t,x,y) + \frac{\Delta x^2}{2}\frac{\partial^2 u}{\partial x^2}(t,x,y) + \frac{\Delta x^3}{6}\frac{\partial^3 u}{\partial x^3}(t,x,y) + \mathcal{O}({\Delta x^4}) $$

ainsi que

$$ u(t,x-\Delta x,y) = u(t,x,y) - \Delta x\frac{\partial u}{\partial x}(t,x,y) + \frac{\Delta x^2}{2}\frac{\partial^2 u}{\partial x^2}(t,x,y) - \frac{\Delta x^3}{6}\frac{\partial^3 u}{\partial x^3}(t,x,y) + \mathcal{O}({\Delta x^4}) $$

En additionnant on obtient

$$ u(t,x-\Delta x,y) + u(t,x+\Delta x,y) = 2u(t,x,y) + \Delta x^2\frac{\partial^2 u}{\partial x^2}(t,x,y) + \mathcal{O}({\Delta x^4}) $$

c'est-à-dire

$$ \frac{\partial^2 u}{\partial x^2}(t,x,y) = \frac{u(t,x-\Delta x,y) -2u(t,x,y) + u(t,x+\Delta x,y)} {\Delta x^2} + \mathcal{O}({\Delta x^2}) $$

On peut alors estimer que

$$\frac{\partial u^2}{\partial x^2}(t_k,x_i,y_j)\approx\frac{U_{i-1,j}^k - 2U_{i,j}^k + U_{i+1,j}^k}{\Delta h^2}$$

On obtient de la même façon l'estimation

$$\frac{\partial u^2}{\partial y^2}(t_k,x_i,y_j)\approx\frac{U_{i,j-1}^k - 2U_{i,j}^k + U_{i,j+1}^k}{\Delta h^2}$$

Ceci conduit à proposer le système d'équations suivant pour les valeurs $U_{i,j}^k$:

$$\frac{U_{i,j}^{k+1}-U_{i,j}^k}{\Delta t} =\alpha \frac{U_{i-1,j}^k - 2U_{i,j}^k + U_{i+1,j}^k + U_{i,j-1}^k - 2U_{i,j}^k + U_{i,j+1}^k}{\Delta h^2} $$

ce qui peut s'écrire, en posant $\lambda = \alpha\Delta t/\Delta h^2$

$$U_{i,\,j}^{k+1} = \lambda\left(U_{i-1,j}^k + U_{i+1,j}^k + U_{i,j-1}^k + U_{i,j+1}^k\right) + (1-4\lambda)U_{i,j}^k \\ \text{pour}\quad k\geq 0\quad\text{et} \quad 1\leq i\leq n-1,\quad 1\leq j\leq n-1 $$

Ce système d'équations doit être complété par les conditions initiales et les conditions aux bords de la plaque :

$$ U_{i,\,j}^0=u_{\text{init}}(x_i,\,y_j) \\ U_{0,\,j}^k=u_{\text{left}}(y_j) \\ U_{n,\,j}^k=u_{\text{right}}(y_j) \\ U_{i,\,0}^k=u_{\text{bottom}}(x_i) \\ U_{i,\,n}^k=u_{\text{top}}(x_i) $$

Consistance et stabilité du modèle discret

Le schéma numérique que nous venons d'établir possède une unique solution mathématique. Une question importante se pose maintenant. Dans quelle mesure la solution exacte de ce schéma numérique est-elle en adéquation avec la solution exacte de l'EDP ? Il serait en particulier souhaitable que la solution exacte de l'EDP aux points du maillage soit une solution approchée du schéma numérique, avec une approximation qui tende vers 0 lorsque les pas $\Delta h$ et $\Delta t$ tendent vers 0. On peut montrer que le schéma numérique que nous avons défini vérifie bien cette propriété. On dit que le schéma numérique est consistant à l'EDP.

Une autre difficulté apparait. Nous ne sommes pas en mesure de déterminer la solution exacte du schéma numérique. Les calculs, effectués en nombres flottants, seront nécessairement approximatifs. De plus, lors du processus itératif, les erreurs d'arrondis peuvent éventuellement s'accumuler et conduire alors à des erreurs importantes à partir d'un certain nombre d'itérations. On dit que le schéma numérique est stable lorsque les perturbations de la solution numérique ne sont pas amplifiées au cours des itérations. On peut montrer que le schéma numérique ci-dessus est stable à condition que $\lambda$ soit inférieur à une certaine valeur que vous aurez à déterminer expérimentalement.

Codage de la fonction simul

La fonction ci-dessous permet d'afficher les paramètres de la simulation. Elle doit être appeler par la fonction de simul ci-dessous pour pouvoir facilement visualiser les paramètres utilisés.

In [ ]:
def print_params(α, h, tmax, n, Δh, s, Δt, λ):
    """ Affiche les paramètres de la simulation
        de façon agréable à lire.
    """
    print("SIMULATION :")
    print("\tDonnées physiques :")
    print("\t\t%-30s : α = %.3g" % ("Coefficient de diffusion", α))
    print("\t\t%-30s : h = %.3g" % ("Largeur de la plaque", h))
    print("\t\t%-30s : tmax = %.3g" % ("Durée de simulation", tmax))
    print("\tDiscrétisation :")
    print("\t\t%-30s : n = %d" % ("Nombre de points", n))
    print("\t\t%-30s : Δh = %.3g" % ("Écart inter-points", Δh))
    print("\t\t%-30s : s = %d" % ("Nombre d'instants", s))
    print("\t\t%-30s : Δt = %.3g" 
            % ("Intervalle entre deux instants", Δt))
    print("\t\t%-30s : λ = %.3g" % ("Coefficient λ", λ))
    print("\tNombre d'opérations élémentaires"
                 ": %.3g millions"
                  % (s * n**2 / 10**6))
    print()

Question 1

Écrire le début de la fonction simul.

La largeur de la plaque est prise pour unité de longueur. Autrement dit $h=1$.

Par exemple

simul(α=0.001, tmax=100, n=80, λ=0.2)

affiche :

SIMULATION :
    Données physiques :
        Coefficient de diffusion       : α = 0.001
        Largeur de la plaque           : h = 1
        Durée de simulation            : tmax = 100
    Discrétisation :
        Nombre de points               : n = 80
        Écart inter-points             : Δh = 0.0125
        Nombre d'instants              : s = 3199
        Intervalle entre deux instants : Δt = 0.0313
        Coefficient λ                  : λ = 0.2
    Nombre d'opérations élémentaires: 20.5 millions

et

simul(α=0.001, tmax=100, n=80, λ=0.8)

affiche

SIMULATION :
    Données physiques :
        Coefficient de diffusion       : α = 0.001
        Largeur de la plaque           : h = 1
        Durée de simulation            : tmax = 100
    Discrétisation :
        Nombre de points               : n = 80
        Écart inter-points             : Δh = 0.0125
        Nombre d'instants              : s = 799
        Intervalle entre deux instants : Δt = 0.125
        Coefficient λ                  : λ = 0.8
    Nombre d'opérations élémentaires: 5.11 millions
In [ ]:
def simul(α, tmax, n, λ):
    """ Calcule et affiche tous les paramètres de la simulation :
        α, h, tmax, n, Δh, s, Δt, λ
    """
    ### À COMPLÉTER
    print_params(α, h, tmax, n, Δh, s, Δt, λ)

Question 2

Compléter le code de la fonction simul en créant et initialisant un tableau Numpy U[k, i, j] avec les conditions initiales et les conditions sur les bords.

Par exemple

U = simul(α=0.001, tmax=100, n=5, λ=0.2,
          u_init=lambda x, y : x + y,
          u_left=lambda y: y,
          u_right=lambda y: 1 - y,
          u_bottom=lambda x: x,
          u_top=lambda x: 1 - x)

print(U[0])

affiche :

SIMULATION :
    Données physiques :
        Coefficient de diffusion       : α = 0.001
        Largeur de la plaque           : h = 1
        Durée de simulation            : tmax = 100
    Discrétisation :
        Nombre de points               : n = 5
        Écart inter-points             : Δh = 0.2
        Nombre d'instants              : s = 12
        Intervalle entre deux instants : Δt = 8
        Coefficient λ                  : λ = 0.2
    Nombre d'opérations élémentaires: 0.0003 millions

[[ 0.   0.2  0.4  0.6  0.8  1. ]
 [ 0.2  0.4  0.6  0.8  1.   0.8]
 [ 0.4  0.6  0.8  1.   1.2  0.6]
 [ 0.6  0.8  1.   1.2  1.4  0.4]
 [ 0.8  1.   1.2  1.4  1.6  0.2]
 [ 1.   0.8  0.6  0.4  0.2  0. ]]
In [ ]:
def simul(α, tmax, n, λ,
          u_init=lambda x, y : 0,
          u_left=lambda y : 0,
          u_right=lambda y : 0,
          u_top=lambda x : 0,
          u_bottom=lambda x : 0):
    """ Retourne un tableau Numpy U[k, i, j] correctement
        initialisé à partir des conditions à l'instant 0
        et des conditions aux quatre bords de la plaque.
    """
    ### À COMPLÉTER

Question 3

Compléter le code de la fonction simul de sorte que la fonction simule l'évolution de température dans la plaque et retourne le tableau U contenant la solution.

Par exemple :

U = simul(α=0.001, tmax=100, n=5, λ=0.2,
          u_left=lambda y: y,
          u_right=lambda y: 1 - y,
          u_bottom=lambda x: x,
          u_top=lambda x: 1 - x)

print(U[2])

affiche :

SIMULATION :
    Données physiques :
        Coefficient de diffusion       : α = 0.001
        Largeur de la plaque           : h = 1
        Durée de simulation            : tmax = 100
    Discrétisation :
        Nombre de points               : n = 5
        Écart inter-points             : Δh = 0.2
        Nombre d'instants              : s = 12
        Intervalle entre deux instants : Δt = 8
        Coefficient λ                  : λ = 0.2
    Nombre d'opérations élémentaires: 0.0003 millions

[[ 0.     0.2    0.4    0.6    0.8    1.   ]
 [ 0.2    0.128  0.136  0.224  0.432  0.8  ]
 [ 0.4    0.136  0.032  0.048  0.224  0.6  ]
 [ 0.6    0.224  0.048  0.032  0.136  0.4  ]
 [ 0.8    0.432  0.224  0.136  0.128  0.2  ]
 [ 1.     0.8    0.6    0.4    0.2    0.   ]]
In [ ]:
def simul(α, tmax, n, λ,
          u_init=lambda x, y : 0,
          u_left=lambda y : 0,
          u_right=lambda y : 0,
          u_top=lambda x : 0,
          u_bottom=lambda x : 0):
    """ Simule l'évolution de la température dans la plaque.
        Retourne le tableau U contenant la solution.
    """
    ### À COMPLÉTER

Question 4

Réaliser une simulation où la température initiale est nulle et la température sur les bords est maintenue à une valeur de 10 et vérifier que la température tend bien vers 10 sur toute la plaque.

Choisir $\alpha=0.001$, $n=5$ et $\lambda=0.2$.

In [ ]:
 

Quelques simulations

Question 5

Réaliser une simulation sur une plaque où la température initiale est nulle et la température sur les bords est maintenue à une valeur de 10. Représenter la distribution de température sur la plaque à l'instant $t_{\text{max}}=100$.

Choisir $n=80$ et $\lambda=0.2$ et simuler pour plusieurs valeurs pour le coefficient de difusion $\alpha$.

Utiliser plt.matshow(M) pour représenter graphiquement les valeurs d'un tableau à deux dimensions M selon une échelle de couleur.

In [ ]:
 

Question 6

Réaliser une simulation sur une plaque où la température initiale est nulle et la température sur les bords est maintenue à une valeur de 10 et représenter l'évolution de la température au centre de la plaque.

Choisir $n=50$ et $\lambda=0.2$ et simuler pour plusieurs valeurs du coefficient de difusion.

In [ ]:
 

Question 7

Réaliser une simulation pour une plaque de coefficient $\alpha=0.001$ où la température initiale est nulle et la température sur les bords est maintenue à une valeur de 10 et représenter l'évolution de la température sur l'axe de symétrie horizontal de la plaque. Pour cela tracer sur le même graphique plusieurs courbes pour différents instants $t$.

Choisir $n=50$ et $\lambda=0.2$.

In [ ]:
 

Stabilité du schéma numérique

Question 8

Dans toutes les simulations précédentes nous avons donné au coefficient $\lambda$ la valeur $0.2$. Ce coefficient correspond au coefficient de proportiion entre $\Delta t$ et $\Delta h^2$. Intuitivement on conçoit bien que le pas de temps $\Delta t$ et le pas d'espace $\Delta h$ doivent être cohérents. En particulier, si $\Delta t$ est trop grand par rapport à $\Delta h$ alors les calculs sont incohérents et le schéma numérique devient instable.

Déterminer expérimentalement la valeur limite de $\lambda$ au delà de laquelle le schéma numérique devient instable.

Expérimenter avec une température initiale nulle, en maintenant le bord gauche à 10 et les autres bords à 0.

Choisir $\alpha=0.001$, $t_{\text{max}}=100$ et $n=80$.

Observer la température finale sur l'axe de symétrie horizontale de la plaque.

In [ ]:
 

Conditions de Von Neumann

Question 9

Les conditions de von Neumann consistent à imposer la valeur de la dérivée de $u$ sur les bords de la plaque. On note ces dérivées $v_{\text{left}}(y)$, $v_{\text{right}}(y)$, $v_{\text{bottom}}(x)$ et $v_{\text{top}}(x)$.

Traduire les conditions de von Neumann par des équations que doivent satisfaire les $U_{i,\,j}^k$.

Question 10

La prise en compte des conditions de Neumann consiste à calculer à la fin de chaque itération les nouvelles valeurs de $U_{0,\,j}^{k+1}$, $U_{n,\,j}^{k+1}$, $U_{i,\,0}^{k+1}$ et $U_{i,\,n}^{k+1}$.

Modifier la fonction simul de sorte qu'elle permette de simuler des conditions de von Neumann.

Par exemple

def u_init(x, y):
    if x < 0.5:
        u1 = x
    else:
        u1 = 1 - x
    if y < 0.5:
        u2 = y
    else:
        u2 = 1 - y
    return 10 * u1 * u2

U = simul(α=0.001, tmax=100, n=5, λ=0.2,
          u_init=u_init,
          v_left=lambda y: 1,
          v_right=lambda y: 1,
          v_bottom=lambda x: 1,
          v_top=lambda x: 1)

print(U[2])

affiche :

SIMULATION :
    Données physiques :
        Coefficient de diffusion       : α = 0.001
        Largeur de la plaque           : h = 1
        Durée de simulation            : tmax = 100
    Discrétisation :
        Nombre de points               : n = 5
        Écart inter-points             : Δh = 0.2
        Nombre d'instants              : s = 12
        Intervalle entre deux instants : Δt = 8
        Coefficient λ                  : λ = 0.2
    Nombre d'opérations élémentaires: 0.0003 millions

[[ 0.048  0.248  0.528  0.528  0.248  0.048]
 [ 0.248  0.448  0.728  0.728  0.448  0.248]
 [ 0.528  0.728  1.056  1.056  0.728  0.528]
 [ 0.528  0.728  1.056  1.056  0.728  0.528]
 [ 0.328  0.528  0.808  0.808  0.528  0.328]
 [ 0.528  0.728  1.008  1.008  0.728  0.528]]
In [ ]:
def simul(α, tmax, n, λ,
          u_init=lambda x, y : 0,
          u_left=lambda y : 0,
          u_right=lambda y : 0,
          u_top=lambda x : 0,
          u_bottom=lambda x : 0,
          v_left=None,
          v_right=None,
          v_top=None,
          v_bottom=None):
    """ Simule l'évolution de la température dans la plaque.
        Retourne le tableau U contenant la solution.
        
        Par chaque bord individuellement :
         * si v_bord vaut None alors la condition
           u_bord est maintenu à tout instant t
         * sinon v_bord est la fonction donnant la
           dérivée à maintenir sur ce bord.
    """    
    ### À COMPLÉTER

Question 11

La température initiale de la plaque de coefficient $\alpha=0.001$ est donnée par la fonction

$$u_{\text{init}}(x,y)= 10\,sin(\pi x)\, sin(\pi y)$$

Les bords de la plaque sont calorifugées. Quelles conditions faut-il imposer aux bords pour modéliser le problème ?

Question 12

Simuler le problème et représenter la distribution de la température à différents instants.

In [ ]: