Différences finies 1D

L'objectif ici est de présenter la discrétisation différence finie 1D pour le problème de Dirichlet homogène: $$ \begin{align} &-u''(x)=f(x),\quad \forall x\in ]0,1[,\\ &u(0)=u(1)=0. \end {align} $$

Approximation de la dérivée seconde

On rappelle que si $u$ est une fonction $\mathcal C^4$, alors $$ -u''(x)=\frac {-u(x-h)+2u(x)-u(x+h)}{h^2}+\mathcal O(h^2). $$ Cette formule, basée sur des évaluations de $u$ pour approcher $u''$ lorsque le réel $h>0$ est petit, va servir d'approximation évaluée en de nombreux points distants de $h$. Il s'en suivra un système linéaire.

On introduit la discrétisation uniforme du segment $]0,1[$ par $$ x_i=i h, \quad i=0\cdots N+1, $$ avec $h=\frac 1 {N+1}$. Pour $i=1$ à $N$, on va approcher l'équation $-u''(x_i)=f(x_i)$. Pour $i=0$ et $i=N+1$, on va utiliser la condition limite. On note $u_i\sim u(x_i)$ et $f_i=f(x_i)$ satisfaisant $$ \begin {align} &u_0=0\\ &\frac {-u_{0}+2u_1-u_2}{h^2}=f_1\\ &\frac {-u_{1}+2u_2-u_3}{h^2}=f_2\\ &\cdots\\ &\frac {-u_{i-1}+2u_i-u_{i+1}}{h^2}=f_i\\ &\cdots\\ &\frac {-u_{N-1}+2u_N-u_{N+1}}{h^2}=f_ N\\ &u_{N+1}=0. \end {align} $$ En éliminant la première et la dernière equation substituées dans la deuxième et avant-dernière equation, il vient le système de $N$ inconnues $U=(u_1,\cdots,u_N)^t$ suivant: $$ A U=F,\quad F=(f_1,\cdots,f_N)^t, $$ avec $$ A=\frac 1 {h^2}\left (\matrix{&2&-1&0&\cdots&0\\ &-1&2&-1&0&\cdots\\ &0&-1&2&-1&0 \\ &-&-&-&-&-\\ &0&\cdots&0&-1&2}\right ) $$ L'objectif est alors de construire ce système linéaire en python. On utilisera la gestion creuse des matrices à l'aide des outils de la bibliothèque scipy.

Dirichlet inhomogène

Construire la solution du problème de Dirichlet inhomogène $$ \begin{align} &-u''(x)=f(x),\quad \forall x\in ]0,1[,\\ &u(0)=a\\ &u(1)=b. \end {align} $$ Choisir $f$ $a$ et $b$ à votre guise.

Condition limite périodique

Soit le problème, $$ \begin{align} &-u''(x)+u(x)=f(x),\quad \forall x\in ]0,1[,\\ &u(0)=u(1), \quad u'(0)=u'(1). \end {align} $$ qui revient à supposer $u$ $1-$périodique et $\mathcal C^1$.

Neuman homogène

Construisons la solution du problème de Neuman homogène $$ \begin{align} &-u''(x) +u(x)=f(x),\quad \forall x\in ]0,1[,\\ &u'(0)=u'(1)=0. \end {align} $$ On introduit la discrétisation uniforme du segment $]0,1[$ par $$ x_i=i h, \quad i=0\cdots N-1, $$ avec $h=\frac 1 {N-1}$. Pour $i=0$ à $N-1$, on va approcher l'équation $-u''(x_i)+u(x_i)=f(x_i)$ de la même manière que pour le problème de Dirichlet. On note $u_i\sim u(x_i)$ et $f_i=f(x_i)$. Pour $i=0$ et $i=N-1$, on introduit alors $u_{-1}$ et $u_{N}$ qu'il convient d'éliminer par la condition limite de Neuman homogène. On remarque que la différence centrée qui approche la dérivée ($u'(x)=\frac {u(x+h)-u(x-h)} {2h}$) permet d'écrire $$ u_{-1}=u_1,\quad u_{N}=u_{N-2}. $$ On a alors le système linéaire $$ \left \{\begin {align} &u_{-1}=u_1\\ &\frac {-u_{-1}+2u_0-u_1}{h^2}+u_0=f_0\\ &\frac {-u_{0}+2u_1-u_2}{h^2}+u_1=f_1\\ &\frac {-u_{1}+2u_2-u_3}{h^2}+u_2=f_2\\ &\cdots\\ &\frac {-u_{i-1}+2u_i-u_{i+1}}{h^2}+u_i=f_i\\ &\cdots\\ &\frac {-u_{N-2}+2u_{N-1}-u_{N}}{h^2}+u_{N-1}=f_{N-1}\\ &u_{N}=u_{N-2}. \end {align}\right . $$ En éliminant la première et la dernière equation substituées dans la deuxième et avant-dernière equation, il vient le système de $N$ inconnues $U=(u_0,\cdots,u_{N-1})^t$ suivant: $$ A U=F,\quad F=\left (f_0,\cdots,f_{N-1}\right)^t, $$ avec $$ A=\frac 1 {h^2}\left (\matrix{&2&-2&0&\cdots&0\\ &-1&2&-1&0&\cdots\\ &0&-1&2&-1&0 \\ &-&-&-&-&-\\ &0&\cdots&0&-2&2}\right )+Id $$ Par soucis de symétrie dans le système et parce que cela est pertinent pour les méthodes de résolution de système linéaire, on symétrise le système en divisant par $2$ la première et dernière équation du système. On obtient alors le système: $$ A U=F,\quad F=\left(\frac {f_{0}}2,f_1,\cdots,f_{N-2},\frac {f_{N-1}}2\right)^t, $$ avec $$ A=\frac 1 {h^2}\left (\matrix{&1&-1&0&\cdots&0\\ &-1&2&-1&0&\cdots\\ &0&-1&2&-1&0 \\ &-&-&-&-&-\\ &0&\cdots&0&-1&1}\right )+\left (\matrix{&\frac 1 2&0& &\cdots&0\\ &0&1&0& &\cdots\\ &&0&1&0& \\ &-&-&-&-&-\\ & &\cdots& &0&\frac 1 2}\right ) $$

L'objectif est alors de construire ce système linéaire en python. On utilisera la gestion creuse des matrices à l'aide des outils de la bibliothèque scipy.

Neuman inhomogène

Construire la solution du problème de Neuman inhomogène $$ \begin{align} &-u''(x) +u(x)=f(x),\quad \forall x\in ]0,1[,\\ &u'(0)=a\\ &u'(1)=b. \end {align} $$ Choisir $f$ $a$ et $b$ à votre guise.

Robin inhomogène

Construire la solution du problème de Robin inhomogène $$ \begin{align} &-u''(x) =f(x),\quad \forall x\in ]0,1[,\\ &-u'(0)+d_0 u(0)=a\\ &u'(1)+d_1 u(1)=b. \end {align} $$ Choisir $f$ $d_0>0$ et $d_1>0$ $a$ et $b$ à votre guise.

Construisons la solution du problème de Robin inhomogène en se basant sur la même approche que pour le problème de Neuman. $$ \begin{align} &-u''(x) +u(x)=f(x),\quad \forall x\in ]0,1[,\\ &-u'(0)+d_0 u(0)=a\\ &u'(1)+d_1 u(1)=b. \end {align} $$ On introduit la discrétisation uniforme du segment $]0,1[$ par $$ x_i=i h, \quad i=0\cdots N-1, $$ avec $h=\frac 1 {N-1}$. Pour $i=0$ à $N-1$, on va approcher l'équation $-u''(x_i)+u(x_i)=f(x_i)$. On note $u_i\sim u(x_i)$ et $f_i=f(x_i)$. Pour $i=0$ et $i=N-1$, on introduit alors $u_{-1}$ et $u_{N}$ qu'il convient d'éliminer par la condition limite. On remarque que la différence centrée qui approche la dérivée ($u'(x)=\frac {u(x+h)-u(x-h)} {2h}$) permet d'écrire $$ u_{-1}=u_1 - 2hd_0 u_0 +2ha,\quad u_{N}=u_{N-2}- 2hd_1u_{N-1} +2hb. $$ On a alors le système linéaire $$ \left \{\begin {align} &u_{-1}=u_1- 2hd_0u_0 +2ha\\ &\frac {-u_{-1}+2u_0-u_1}{h^2}+u_0=f_0\\ &\frac {-u_{0}+2u_1-u_2}{h^2}+u_1=f_1\\ &\frac {-u_{1}+2u_2-u_3}{h^2}+u_2=f_2\\ &\cdots\\ &\frac {-u_{i-1}+2u_i-u_{i+1}}{h^2}+u_i=f_i\\ &\cdots\\ &\frac {-u_{N-2}+2u_{N-1}-u_{N}}{h^2}+u_{N-1}=f_{N-1}\\ &u_{N}=u_{N-2}- 2hd_1u_{N-1} +2hb. \end {align}\right . $$ En éliminant la première et la dernière equation substituées dans la deuxième et avant-dernière equation, il vient le système de $N$ inconnues $U=(u_0,\cdots,u_{N-1})^t$ suivant: $$ A U=F,\quad F=\left (f_0,\cdots,f_{N-1}\right)^t+\left (\frac {2a} h,0,\cdots,0,\frac {2b} h\right)^t, $$ avec $$ A=\frac 1 {h^2}\left (\matrix{&2+2hd_0&-2&0&\cdots&0\\ &-1&2&-1&0&\cdots\\ &0&-1&2&-1&0 \\ &-&-&-&-&-\\ &0&\cdots&0&-2&2+2hd_1}\right )+Id $$ Par soucis de symétrie dans le système et parce que cela est pertinent pour les méthodes de résolution de système linéaire, on symétrise le système en divisant par $2$ la première et dernière équation du système. On obtient alors le système: $$ A U=F,\quad F=\left(\frac {f_{0}}2,f_1,\cdots,f_{N-2},\frac {f_{N-1}}2\right)^t+\left (\frac {a} h,0,\cdots,0,\frac {b} h\right)^t, $$ avec $$ A=\frac 1 {h^2}\left (\matrix{&1+hd_0&-1&0&\cdots&0\\ &-1&2&-1&0&\cdots\\ &0&-1&2&-1&0 \\ &-&-&-&-&-\\ &0&\cdots&0&-1&1+hd_1}\right )+\left (\matrix{&\frac 1 2&0& &\cdots&0\\ &0&1&0& &\cdots\\ &&0&1&0& \\ &-&-&-&-&-\\ & &\cdots& &0&\frac 1 2}\right ) $$

L'objectif est alors de construire ce système linéaire en python. On utilisera la gestion creuse des matrices à l'aide des outils de la bibliothèque scipy. On remarque que ce problème contient en particulier le cas de Neuman inhomogène en choisissant $d_0=0$, $d_1=0$.

Estimation d'erreur

On va estimer l'erreur entre la solution approchée par la méthode des différences finies centrées et la solution exacte du problème . La théorie dit que l'erreur d'approximation est d'ordre 2, c'est à dire de taille $\mathcal O(h^2)$ et fait plus précisément intervenir la dérivée quatrième de la solution exacte. Ainsi si la solution est un polynôme de degré inférieur ou égal à 3, la solution approchée coincide avec la solution exacte en chacun des points de la discrétisation $x_i$. Il conviendra donc de choisir une solution exacte qui n'est pas un olynôme de degré inférieur à 3 et qui satisfait bien-sûr les conditions limites.

Problème de Dirichlet

Choisissons par exemple la solution non triviale $u(x)=\sin(x^2)\sin(x-1)$ du problème de Dirichlet $$ \begin{align} &-u''(x)=f(x),\quad \forall x\in ]0,1[,\\ &u(0)=u(1)=0. \end {align} $$ Il convient de calculer $f$ afin de satisfaire que $u(x)=\sin(x^2)\sin(x-1)$ est bien solution. Pour des problèmes plus compliqué, on pourra faire appel au calcul formel avec la librairie sympy.

On sait désormais l'expression de $f$. On va alors construire cette fonction et définir la solution exacte. En définissant l'erreur entre la solution exacte et approchée, pour la norme de notre choix, comme une fonction python qui dépend de la taille de la discrétisation, on pourra alors calculer aisément l'erreur pour plusieurs discrétisation. La fonction construite sera également paramètrée par la fonction $f$, $u$ (solution exacte) et un entier qui caractérise la norme choisie.

Différences finies 2D

L'objectif ici est de présenter la discrétisation différence finie 2D pour le problème de Dirichlet homogène posé sur un rectangle: $$ \begin{align} -\Delta u(x,y)=f(x,y),&\quad \forall (x,y)\in \Omega=]0,L[\times]0,H[,\\ u(x,y)=0,& \quad \forall (x,y)\in \partial \Omega. \end {align} $$

Approximation du Laplacien

On utilise la même approximation que dans le cas 1D pour la direction $x$ et pour la direction $y$, si $u$ est une fonction $\mathcal C^4$, alors $$ -\Delta u(x,y)=\frac {-u(x-l,y)+2u(x,y)-u(x+l,y)}{l^2}+\frac {-u(x,y-h)+2u(x,y)-u(x,y+h)}{h^2}+\mathcal O(l^2+h^2). $$ Cette formule, basée sur des évaluations de $u$ sur une grille Cartésienne de pas $l$ dans la direction $x$ et de pas $h$ dans la direction $y$, va servir d'approximation de l'équation. Il s'en suivra un système linéaire dont le nombre d'inconnues est le nombre de points de grille internes au domaine $\Omega$. On représente les points de la grille pour lesquels l'équation sera discrétisée.

Numérotation

Il convient alors de numéroter ces points à l'aide d'un seul indice en choisissant un sens de parcours de tous ces points. Puis, on écrira pour chacun d'eux la discrétisation et ainsi les coefficients de la matrice correspondante à ce parcours, ie cette numérotation. Pour le point numéro $k$, on sera en mesure de remplir la ligne numéro $k$ de la matrice. Celle-ci contiendra donc $NM$ lignes et colonnes.

La manière la plus classique de numéroter est de parcourir de gauche à droite depuis la ligne du bas jusqu'à la ligne du haut. Le premier point de chaque ligne est un point particulier car son voisin gauche tombe sur le bord. Il en va de même du dernier de chaque ligne dont le voisin de droite est sur le bord droit. La première et dernière ligne sont aussi des cas particuliers puisque le voisin du bas ou du haut tombe sur la bord bas ou haut.

Si on considère le noeud $k$ associé aux indices $1\le i\le N$ et $1\le j\le M$ de la numérotation horizontale et verticale, on a $k=num(i,j)=i+(j-1)N$. Si on choisit $k=num(i,j)$ de sorte que $1< i< N$ et $1< j< M$ (noeud non voisin du bord), le voisin de gauche a le numéro $k-1$, le voisin de droite a le numéro $k+1$, le voisin du bas a le numéro $k-N$, le voisin du haut a le numéro $k+N$.

Structure de la matrice

On en déduit que la matrice sera pentadiagonale avec une diagonale principale, une sous-diagonale décalée de $1$, une sous-diagonale décalée de $N$, une sur-diagonale décalée de $1$, une sur-diagonale décalée de $N$.
Il faudra gérer les cas particuliers des noeuds proche du bord afin de prendre en compte les voisins du bord pour lequel la solution est nulle (Dirichlet inhomogène).

L'assemblage se construit donc à l'aide de 5 diagonales dont il faudra corriger ou vérifier les $N$ premières lignes, les $N$ dernières lignes, ainsi que la première ligne de chaque sous-bloc de taille $N$ (points gauches) et ainsi que la dernière ligne de chaque sous-bloc de taille $N$ (points droits).

On doit se convaincre sur un brouillon que la matrice obtenue est symétrique.

Remplissage de la matrice en python

Validation Dirichlet inhomogène

Coder le problème de Dirichlet inhomogène et le valider On construit une validation du problème de Dirichlet inhomogène.

On choisit une fonction non triviale , y compris sur le bord. On calcule le second membre correspondant. On reprend la construction de la solution approchée de Dirichlet homogène et on vient rajouter dans un vecteur au second membre la contribution des termes de Dirichlet inhomogène. On évalue alors la solution exacte en ces points du bord puisque ce sont une donnée du problème.

Puis on calcule l'erreur entre la solution exacte et approchée. Pour se simplifier la vie, on lie $N$ et $M$ pour regarder l'erreur comme fonction d'un paramètre, $N$.

Robin inhomogène

Coder le problème de Robin inhomogène à coefficients variables, puis le valider. Il est conseillé de procéder par étape, en construisant et validant d'abord Robin homogène à coefficient constant.

Après avoir compris et valider un problème de Robin 1D, l'extension 2D est naturelle, à condition d'avoir compris la complexité de la numérotation en 2D. La grille de discrétisation couvre les bords du domaine. Sur les noeuds du bord, la condition limite est utilisée pour éliminer le point fictif intervenant dans le schéma à 5 points. Comme en 1D, une approximation centrée de la dérivée permet de conserver des approximations à erreur quadratique.

Les 4 points de la discrétisation correspondant aux 4 coins font intervenir 2 points fictifs, ils sont éliminés en utilisant les CL pour la normale horizontale (terme en $\partial_x^2$) et pour la normale verticale (terme en $\partial_y^2$). On utilise la normale qui nous arange selon le terme vu que la normale n'est pas définie en ce point.