On s'intéresse à la résolution numérique du problème de valeurs initiales et aux limites pour l'équation de la chaleur, avec conditions aux limites de Dirichlet, sur un domaine $]a,b[\times[0,T]$ :
\begin{equation} (C)\ \ \ \ \ \begin{cases} \displaystyle{\frac{\partial u}{\partial t}(t,x)-\frac{\partial^2 u} {\partial x^2}(t,x)=f(t,x)},&x\in]a,b[,\ t\le T,\\[4pt] u(0,x)=u^0(x),& x\in]a,b[,\\[4pt] u(t,a)=g_a(t),\ u(t,b)=g_b(t),&t\le T, \end{cases} \end{equation}
où $u^0:[a,b]\longrightarrow\RR$ est une fonction de classe $C^2$ donnée, et $g_a,\ g_b:[0,+\infty[\longrightarrow\RR$ sont des fonctions de classe $C^1$ données, vérifiant $g_a(0)=u_0(a),\ g_b(0)=u_0(b)$.
Dans cette équation l'inconnue $u$ représente par exemple la distribution de chaleur à l'instant $t$ en chaque point $x$ d'un matériau homogène conducteur de la chaleur, que l'on assimile à un domaine $\Omega$ de $\RR^n,\ n=1,\ 2$ ou $3$ (ici pour simplifier on se restreindra à la dimension $n=1$). La fonction $f$ modélise une source de chaleur extérieure à laquelle est soumis le matériau.
On suppose en plus que la température est fixe au bord, ce qui se traduit par des conditions aux limites de Dirichlet, on aurait pu aussi supposer par exemple que le flux de chaleur au bord dans la direction normale est nul, ce qui se traduirait par une condition de Neumann.
Notre objectif est de calculer des valeurs approchées de la solution exacte $u$ de $(C)$ par la méthode des différences finies. Pour cela on introduit :
pour $M>0$ donné, une discrétisation uniforme de l'intervalle $[a,b]$ définie par les $M+2$ points $x_j=a+jh,\ j=0,\dots,M+1$, avec $h=\frac{b-a}{M+1}$ le pas de la discrétisation (aussi appelé pas d'espace) ;
le pas de temps $k>0$ et les instants temporels $t^n=nk,\ n\in\N,n\leq N,$ avec $Nk=T$.
On cherche alors des valeurs $u^n_j$ approchant $u(t^n,x_j)$, pour $n\in\N,n\leq N$ et $j=0,\dots,M+1$. Comme la solution exacte $u$ de $(C)$ vérifie
$$ u(0,x_j)=u_0(x_j), $$
il est naturel de prendre $u^0_j=u_0(x_j),\ j=0,\dots,M+1$.
On considèrera deux schémas pour calculer les valeurs approchées $(u^n_j)_{j=0,\dots,M+1}\,,\ \ $ pour $n=1,\dots,N$.
Le schéma explicite : \begin{equation*} \begin{cases} \displaystyle{\frac{u^{n+1}_j-u^n_j}{k}-\frac{u_{j+1}^{n}-2u_j^{n}+u_{j-1}^{n}}{h^2}=f(t^n,x_j)},&j=1,\dots,M,\ n=0,\dots,N-1,\\[4pt] u^{n+1}_0=g_a(t^{n+1}),\ u^{n+1}_{M+1}=g_b(t^{n+1}),&n=0,\dots,N-1. \end{cases} \end{equation*}
Le schéma implicite : \begin{equation*} \begin{cases} \displaystyle{\frac{u^{n+1}_j-u^n_j}{k}-\frac{u_{j+1}^{n+1}-2u_j^{n+1}+u_{j-1}^{n+1}}{h^2}=f(t^{n+1},x_j)},&j=1,\dots,M,\ n=0,\dots,N-1,\\[4pt] u^{n+1}_0=g_a(t^{n+1}),\ u^{n+1}_{M+1}=g_b(t^{n+1}),&n=0,\dots,N-1. \end{cases} \end{equation*}
Nous allons comparer les résultats des deux schémas précédents dans le cas de l'équation de la chaleur homogène ($f=0$), avec conditions aux limites de Dirichlet homogènes ($g_a=g_b=0$). En particulier nous allons voir comment les deux schémas sont influencés par la valeur du paramètre
$$ \lambda=\frac{k}{h^2}. $$
Question 1.
Écrire une fonction python DFchaleurExp(a,b,M,T,lambda,u_0) qui calcule et retourne la discrétisation en espace x et la solution approchée U de $(C)$ avec $f=g_a=g_b=0$, à un temps final $T$ pour une valeur $\lambda$ fixée, obtenue par la méthode explicite.
On prendra garde à ne stocker que la solution au temps final dans le
programme. On n'oubliera pas de définir la valeur approchée en $x=a$ et en $x=b,\ u_{0}^{n}$ et $u^n_{M+1}$.
Question 2. En notant $U^n = (u_1^n,\dots,u^n_M)$, montrer que le schéma implicite peut être mis sous la forme matricielle
$$ \left(\mathrm{I}_M - k A_h \right)U^{n+1} = U^n, $$
où la matrice $A_h$ est à déterminer. Écrire une fonction python DFchaleurImp(a,b,M,T,lambda,u_0) qui calcule et retourne la discrétisation en espace x et la solution approchée U de $(C)$ avec $f=g_a=g_b=0$, à un temps final $T$ pour une valeur $\lambda$ fixée, obtenue par la méthode implicite.
Question 3. On considère $u^0(x)=\sin(\pi x)$ et $[a,b]=[0,1].$ Vérifiez que la solution de (C) est dans ce cas $u(t,x)=e^{-{\pi^2 t}}\sin(\pi x).$
Pour $M=49,$ calculer la solution approchée donnée par la méthode explicite à l'instant $T=2,$ et représenter sur un même graphe la solution approchée et la solution exacte, en utilisant $\lambda=\frac14$. Que se passe-t-il si on considère $\lambda>\frac12$ ?
Faire de mếme pour la méthode implicite. Quelles différences constatez-vous ?
Question 4. Adaptez vos programmes pour prendre en compte un second membre non nul, des conditions aux limites non homogènes, des conditions de Neumann, homogènes ou non.
Supposons ici que la solution $u$ de $(C)$ est une fonction de classe $C^2$ en temps et de classe $C^4$ en espace.
Nous souhaitons analyser la convergence des schémas étudiés.
On définit ainsi :
L'erreur locale $e^n_j$ à l'instant $t^n$, au point $x_j$, par $e^n_j=u(t^n,x_j)-u^n_j,\ $ pour $0\leq n\leq N$ et $0\leq j\leq M+1$.
L'erreur de consistance locale $\varepsilon^n_j(u)$ relative à la solution $u$ en $(t^n,x_j)$ par
$$ \varepsilon^n_j(u)=\frac{u(t^{n+1},x_j)-u(t^n,x_{j})}{k}-\frac{u(t^n,x_{j+1})-2u(t^n,x_{j})+u(t^n,x_{j-1})}{h^2}, $$ pour $1\leq j\leq M$.
Nous voudrons montrer que le schéma explicite converge à l'ordre 1 en temps et 2 en espace, sous la condition
$$ \frac{k}{h^2}\leq\frac12, $$
plus précisément nous voudrons montrer que si la condition ci-dessus est vérifiée, il existe $C>0$, indépendante de $k$ et de $h$, tel que
$$ \sup_{n\in\N,\ nk\leq T}\left(\max_{j=0,\dots,M+1}|e^n_j|\right)\leq C(k+h^2) $$
et le schéma converge donc au sens où
$$ \lim_{h,k\to 0,\ \frac{k}{h^2}\leq\frac12}\,\sup_{n\in\N,\ nk\leq T}\left(\max_{j=0,\dots,M+1}|e^n_j|\right)=0. $$
Nous pouvons montrer le même type de résultat pour le schéma implicite, sans contrainte sur les pas de temps et d'espace.
Question 1. Montrer que la suite des erreurs $(e^n_j)_{0\leq j\leq M+1,\ 0\leq n\leq N}$ vérifie $e^0_j=0,\ j=0,\dots,M+1$, et, pour $0\leq n\leq N-1$, $$ \begin{cases} e^{n+1}_j=\displaystyle{e^n_{j}+\frac{k}{h^2}(e^n_{j+1}-2e^n_{j}+e^n_{j-1})+k\varepsilon^n_j(u)},\ j=1,\dots,M,\\[8pt] e^{n+1}_0=e^{n+1}_{M+1}=0. \end{cases} $$
Question 2. Montrer que si la condition $\frac{k}{h^2}\leq\frac12$ est vérifiée, $(e^n)_n$ vérifie $$ \|e^{n+1}\|_{\infty}\leq\|e^{n}\|_{\infty}+k\|\varepsilon^n(u)\|_{\infty},\ n=0,\dots,N-1. $$
La condition $\frac{k}{h^2}\leq\frac12$ est une condition de stabilité. On dit que le schéma est* stable sous la condition $\frac{k}{h^2}\leq\frac12$.
Question 3. Montrer qu'il existe $C>0$, indépendant de $h$ et de $k$, tel que pour tout $n\in{0,\dots,N-1},$ on ait $$ \|\varepsilon^n(u)\|_\infty\leq C(h^2+k). $$
On dit que le schéma est consistant d'ordre 2 en espace et d'ordre 1 en temps.
Question 4. En déduire que pour tout $n\in\{0,\dots,N\}$, $$ \|e^n\|_{\infty}\leq T C(h^2+k). $$
et conclure que $$ \lim_{h,k\to 0,\ \frac{k}{h^2}\leq\frac12}\,\sup_{n\in\N,\ nk\leq T}\left(\max_{j=0,\dots,M+1}|e^n_j|\right)=0. $$
Question 5. Montrer le même résultat pour le schéma implicite, sans condition sur $h$ et $k$, en montrant que pour la suite $(u^n_j)_{n,j}$ définie par le schéma implicite on a $$ \|e^{n+1}\|_{\infty}\leq\|e^{n}\|_{\infty}+k\|\varepsilon^n(u)\|_{\infty},\ n=0,\dots,N-1. $$
Pour ce faire, considérer $j_0$ tel que $u^{n+1}_{j_0}=\max_{j}|u^{n+1}_j$.
Le schéma de Cranck Nicolson. Il s'agit d'un schéma d'ordre 2 en temps et en espace, défini par l'algorithme
\begin{equation*} \displaystyle{\frac{u^{n+1}_j-u^n_j}{k}-\frac12\frac{u_{j+1}^{n}-2u_j^{n}+u_{j-1}^{n}}{h^2}-\frac12\frac{u_{j+1}^{n+1}-2u_j^{n+1}+u_{j-1}^{n+1}}{h^2}=f(t^n,x_j)},\ \ j=1,\dots,M,\ n=0,\dots,N-1, \end{equation*}
On peut vérifier qu'il est en effet d'ordre 2 en temps et en espace.