Résolution de l'équation stationnaire de la chaleur en 1-D

Le but de ce TP est de calculer la répartition de température dans une barre. On se place dans le cas où le temps est considéré comme n'intervenant plus dans le phénomène de diffusion. Le problème devient alors stationnaire. L'équation stationnaire de la chaleur en 1-D sera modélisée à l'aide d'une équation aux dérivées partielles avec les conditions aux limites suivantes:

où f désigne une source de chaleur. Il s'agit donc d'une équation de Poisson dans le cas d'une source non-nulle, et d'une équation de Laplace dans le cas contraire.

Avant de commencer

Question 1 : Formation du système linéaire

Nous commençons par construire un maillage du domaine, c'est-à-dire un recouvrement de Ω par des segments de taille Δx. La réunion des segments est égale à Ω, l'intersection de deux segments distincts est soit vide, soit un sommet commun. On utilise le schéma différences finies suivant:

avec N le nombre de points et le pas d'espace vaut alors

Du fait des conditions de Dirichlet, l'inconnue est discrétisée sur les points intérieurs uniquement en adoptant la convention :

On note les vecteurs U et F:

Le vecteur U est alors solution du système linéaire suivant:

avec la matrice tridiagonale A suivante

Cette matrice est creuse, contient trois éléments non-nuls par ligne (sauf sur la première et dernière ligne). A est symétrique définite positive (ses valeurs propres sont strictement positives).

Question 1.1

Implémenter une sous-routine qui génère la matrice A pour un N donné en argument. Vérifier que la matrice est correcte pour N=2, 3 et 4.

Résolution du système linéaire

On cherche maintenant la solution du système

A est symétrique définie positive, donc elle peut être mise sous la forme du produit d'une matrice triangulaire inférieure L et de sa transposée LT :

Pour ce faire, nous allons utiliser la factorisation de Cholesky (vous pouvez aussi regarder votre cours).

Question 2.1

Implémenter une sous-routine effectuant la décomposition de Cholesky de la matrice A. Tester cette sous-routine en vérifiant que la norme de A - L LT est de l'ordre de la précision machine pour différentes valeurs de N (on pourra utiliser les fonctions Fortran matmul, transpose et maxval). Pour débugguer, on pourra éventuellement comparer avec le L produit avec ipython --pylab:

N = 4
A = 2*(N+1)**2*diag(ones(N)) - (N+1)**2*diag(ones(N-1), -1) - (N+1)**2*diag(ones(N-1), 1)
L = cholesky(A)
L

Question 2.2

On veut résoudre

à partir du facteur L précédemment calculée. On doit donc résoudre deux systèmes triangulaires successifs:

Écrire une sous-routine calculant U, L et F étant donnés en argument de la sous-routine. On pourra ne donner qu'un seul vecteur en argument (qui contiendra F en entrée et U en sortie). On n'a ici pas besoin de stocker un vecteur y puisqu'on peut écraser au fur et à mesure les valeurs de F jusqu'à obtenir la solution U. Tester cette sous-routine avec un second membre aléatoire (on utilisera la fonction rand). Vérifier que la norme de A U - F est de l'ordre de la précision machine. Pour débugguer avec ipython, on pourra écrire dans ipython --pylab:

f = rand(N)
y  = linalg.solve(L, f)
u = linalg.solve(L.T, y)
u

Question 3

Question 3.1

Écrire la solution U dans un fichier texte. Dans ipython --pylab, on chargera la solution en tapant:

u = loadtxt("u.dat")

Visualiser le résultat pour

les N/2 premières valeurs de F sont égales à1. Visualiser le résultat (avec plot) pour N=100 et N=1000.

Optimisation de la résolution numérique

Comme A est une matrice tridiagonale, la matrice L est bidiagonale. En conséquence, il est plus économique de stocker A sous la forme de vecteurs d et c:

Le facteur L sera aussi stocké avec deux vecteurs :

Question 4.1

Reprendre les sous-routines de création de A, factorisation et résolution pour ne prendre en argument que deux vecteurs (au lieu d'une matrice). On adaptera les algorithmes en utilisant la relation suivante:

On pourra créer deux programmes distincts. Dans la version tridiagonale, on ne fera pas de vérification (on vérifiera qu'on obtient les mêmes valeurs que dans la version pleine).

Question 4.2

Mesurer et comparer les temps CPU de la décomposition pleine et tridiagonale de la matrice A. Pour mesurer le temps CPU d'une portion de votre code, vous pouvez utiliser la fonction CPU_TIME de la manière suivante:

real*8 :: t1, t2

call CPU_TIME(t1)

// ici le code dont on veut mesurer le temps

call CPU_TIME(t2)

print*, "Elapsed time = ", t2-t1

Analyse de la qualité des résultats numériques

On cherche souvent à évaluer l'erreur numérique commise par une méthode sous la forme :

On appelle p l'ordre de convergence de la méthode. Notez que sous forme logarithmique, l'erreur est linéire (en négligeant le terme en O(1/Np+1) )

Le but de cette question est de quantifier la qualité de la solution numérique calculée à l'aide de la connaissance d'une solution analytique exacte. Considérez par exemple la solution analytique

On remarque que u est nulle en x=0 et x=1 et que

Question 5.1

Le protocole de test est le suivant

Question 5.2

Le but de cette question est de quantifier la qualité de la solution numérique calculée lorsque l'on ne dispose pas de la connaissance d'une solution analytique exacte. Le protocole de test est le suivant:

On pourra tester sur cette méthodologie sur le cas

puis sur le cas de la question 3 :