Skip to content
Extraits de code Groupes Projets
diffusion.tex 8,51 ko
Newer Older
  • Learn to ignore specific revisions
  • Michel Crucifix's avatar
    Michel Crucifix a validé
    
    \section{Equation de diffusion}
    \begin{frame}
    \frametitle{Définition}
    La forme générale de l'\'equation de diffusion s'écrit:
    \begin{equation}
    a(x, t) \pdxx{u} + b(x, t) \pdx{u} + c(x, t) u = \pdt u + f(x, t), \quad
    \text{pour } \left\{ \begin{array}{l} 0 < x < \ell, \\ 0 < t. \end{array} \right..
    \label{eq:chaleur}
    \end{equation}
    Nécessite des conditions initiales \emph{et} des conditions aux frontières. 
    \end{frame}
    \begin{frame}
    \begin{example}
    L'équation de la chaleur (adimensionnelle) se réduit à $a(x, t) = 1$; $b=c=f=0$. \\
    Si on admet, pour conditions ($\ell = 1$):
    $$u(0, t)= u(1, t)=0 \quad \text{ pour }\ t>0\text{, et} $$
    $$u(x, 0) = g(x), \quad \text{pour } 0 \leq x \leq 1, $$
    alors le système admet pour solution \footnote{truc: écrire $u=A(x)C(t)$}:
    \begin{equation}
    u(x, t) = \sum_{n=1}^{\infty} A_n e^{-\lambda_n^2 t} \sin (\lambda_n x),
    \label{}
    \end{equation}
    $\lambda_n=n\pi$ et $A_n = 2 \int_0^{1} g(x) \sin (\lambda x)\, \d x$
    \end{example}
    \end{frame}
    \begin{frame}
    \begin{exercice}
    \label{eqchal}
    Déterminer les solutions analytiques de l'équation de la chaleur adimensionnelle définie page prédente, pour
    \begin{enumerate}
    \item $g(x) = \sin(2\pi x)$, 
    \item $g(x) = \left\{ \begin{array}{ll} 1\quad & \text{si \ } a \leq x \leq b, \\ 0 & sinon. \end{array} \right.$
    \end{enumerate}
    \end{exercice}
    \begin{block}{Propriétés de l'équation de la chaleur}
    \begin{itemize}
    \item La solution est continuement dérivable à l'intérieur des frontières
    \item Les minima et maxima se trouvent aux frontières, soit latérales, soit en $t=0$.
    \item L'information se propage instantanément (illustration: cas 2. de l'exercice ci-dessus)
    \end{itemize}
    \end{block}
    \end{frame}
    \begin{frame}
    \frametitle{Résolution en 4 étapes}
    \begin{enumerate}
    \item Définir les points de discrétisation
    \item Estimer le système en ces points
    \item Le formuler en différences finies
    \item Éliminer les termes de troncature
    \end{enumerate}
    \end{frame}
    \begin{frame}
    \frametitle{1. Discrétisation}
    \begin{tikzpicture}[xscale=1.3]
    \draw [->] (-0.2, 0) -- (4.8, 0) node [yshift=-1em] {$x$} ; 
    \draw [->] (0, -0.2) -- (0, 4.8) node [xshift=-1.5em] {$t$} ;
    \foreach \i in {1, 2, 3}
     { 
     \draw (\i, 0) node [yshift=-1em] {$x_\i$} -- (\i, 4) ; 
     \draw (0, \i) node [xshift=-1.5em] {$t_\i$} -- (4, \i) ; 
     \foreach \j in {1, 2, 3, 4}
      { \node at (\i, \j)  [shape=circle, draw, color=black] {} ; }
      }
      \draw ( 4, 0) node [yshift=-1em] {$\ell$} -- ( 4, 4) ; 
      \draw ( 0, 4) node [xshift= -1.5em] {$t_4$} -- ( 4, 4) ; 
    
      \foreach \j in {0,1, 2, 3,4}
      { \node at (\j, 0 )  [shape=circle, draw, color=black, fill=black!30] {} ; 
        \node at (0 , \j)  [shape=circle, draw, color=black, fill=black!30] {} ; 
        \node at (4 , \j)  [shape=circle, draw, color=black, fill=black!30] {} ; }
    
    \end{tikzpicture}
    
    Pour un grille régulière, on aura $t_j = j \dt$ et $x_i = i h $. 
    \end{frame}
    
    \begin{frame}
    \frametitle{2. : exprimer le système aux points de grille}
    Nous reprenons l'équation de la chaleur, cette fois non-homogène :
    \begin{equation*}
    \pdxx{u} = \pdt{u} + f(x, t),
    \end{equation*}
    avec $u(0, t) = u(1, t) = 0$; $u(x, 0) = g(x)$, et nous adoptons pour notation: \\
    $u_t := \pdt{u}$, $u_x := \pdx{u}$, $u_{xx} := \pdxx{u}$ etc. 
    
    \bigskip
    
    Le système discrétisé est donc : 
    \begin{equation}
    u_{xx} (x_i, t_j) = u_t (x_i, t_j) + f (x_i, t_j)
    \label{}
    \end{equation}
    
    \end{frame}
    \begin{frame}
    \frametitle{3. Exprimer le système sous forme de différence finie}
    On choisit un schéma centré dans l'espace:
    \begin{align*}
    u_{xx}(x_i, t_j) &= \frac{u(x_{i+1}, t_j) - 2u (x_i, t_j) + u(x_{i-1}, t_j)}{h^2} + \mathcal{O}(h^2)
    \end{align*}
    et Euler avant dans le temps:
    \begin{align*}
    u_t(x_i, t_{j}) &= \frac{u(x_i, t_{j+1}) - u (x_i, t_{j})}{k} + \mathcal{O}(k)
    \end{align*}
    \end{frame}
    \begin{frame}
    \frametitle{4. \ldots et on laisse tomber les termes de troncature.}
    Posons $\lambda = k / h^2$, on a:
    \begin{equation}
    u_{i, j+1} = \lambda u_{i+1, j} + (1-2\lambda) u_{i,j} + \lambda u_{i-1, j} - k f_{i,j},
    \label{}
    \end{equation}
    ce qui peut s'écrire sous la forme 
    \begin{align*}
    \vec{u}_0 &= \vec{g}  \\
    \vec{u}_{j+1} &=  (1-\vec{H})\vec{u}_{j} - k \vec{f}_j, \quad \text{pour\ }j = 1\ldots M-1\\
    \end{align*}
    \begin{columns}
    \column{0.35\textwidth}
    où : 
    \begin{align*}
     \vec{u}_j &= (u_{1, j}, u_{2, j}, \ldots , u_{N, j} ) \\
     \vec{f}_j &= (f_{1, j}, f_{2, j}, \ldots , f_{N, j} )  \\
     \vec{g} &= (g_{1}, g_{2}, \ldots , g_{N} )  \\
    \end{align*}
    \column{0.55\textwidth}
    et $\vec{H}=\left( 
    \begin{array}[h]{cccc}
      2\lambda &  -\lambda  &          &           \\
    
    -\lambda    &  2\lambda & -\lambda  &           \\
    
               &\ddots     & \ddots    & \ddots    \\
    
               &           &-\lambda    &  2\lambda   \\
    
    \end{array}
    \right)$
    \end{columns}
    \end{frame}
    \def\stencil{
     \draw [-] (-1.8, 0) node [left] {$t_j$} -- (1.8, 0); 
     \draw [-] (-1.8, 1) node [left] {$t_{j+1}$} -- (1.8, 1); 
     \draw [-] (-1, -0.5) node [yshift=-1.5em] {$x_{i-1}$} -- (-1, 1.5 ); 
     \draw [-] (0, -0.5) node [yshift=-1.5em] {$x_{i}$} -- (0, 1.5 ); 
     \draw [-] (1, -0.5) node [yshift=-1.5em] {$x_{i+1}$} -- (1, 1.5 ); 
     \node at (0,1)  [shape=circle, draw, color=black] { } ; 
     }
    
    
    \begin{frame}
    Le \emph{stencil} représente les points dont dépend le calcul du $u_{j+1, i}$:
    \begin{tikzpicture} 
     \stencil
     \foreach \i in {-1, 0, 1} {\node at (\i,0)  [shape=circle, draw, fill=black!30] { } ; } 
    \end{tikzpicture}
    
    Le système qui nous occupe est donc \emph{explicite.}
    Ce schéma Euler avant de l'équation de la chaleur admet une solution analytique  si $\vec{f}=0$:
    
    \framebox[\textwidth]
    {
    si $\vec{u}_0=w_0 e^{rx\dot\imath}$, alors
    $\vec{u}_n=(1-q)^n \vec{u}_0$, où $q=4\lambda \sin^2\left( \frac{rh}{2} \right)$.
    }
    
    ($\icmplx:= \sqrt{-1}$)
    \end{frame}
    
    \begin{frame}
    
    \renewcommand\arraystretch{2.4} 
    \begin{tabular}{l|l|l}
    %\setlength\minrowclearance{2.4pt}
    \textbf{schéma} & \textbf{Itération} & \textbf{Solution analytique$^\star$ }\\
    \hline
    E. avant & $\vec{u}_{j+1} =  (1-\vec{H})\vec{u}_{j}  $ & $\vec{u}_n=(1-q)^n \vec{u}_0$  \\
    \hline
    E. arrière & $(1+\vec{H})\vec{u}_{j+1} =  \vec{u}_{j} $ & $\vec{u}_n=(1+q)^{-n} \vec{u}_0 $  \\
    \hline
    CN & $(1+\vec{H}/2)\vec{u}_{j+1} =  (1-\vec{H}/2)\vec{u}_{j} $ & $\vec{u}_n=\left(\frac{1-q/2}{1+q/2}\right)^n \vec{u}_0$  \\
    \end{tabular}
    \renewcommand\arraystretch{1.2} 
    \bigskip
    \begin{tabular}{ll}
    E. avant : & Euler avant \\
    E. arrière : & Euler arrière \\
    CN : & Crank - Nicholson 
    \end{tabular}
    
    $\star: $ solution pour $\vec{u}_0 = w_0 e^{rx\icmplx}$ 
    
    \begin{exercice}
    Lesquelles de ces méthodes préservent la propriété de propagation instantanée ?
    \end{exercice}
    \end{frame}
    
    \begin{frame}
    \frametitle{Stabilités A et L}
    Définissons $k(q)$ tel que $\vec{u}_{i+1} = k(q) \vec{u}_{i} $, $\vec{u}_{0}=w_0 e^{\icmplx rx}$:
    \begin{definition}
    La méthode est \emph{A-stable} si $\lim_{i\rightarrow \infty} \vec{u}_i = 0$ \\
    La méthode est \emph{L-stable} si $\lim_{q\rightarrow \infty} k(q) = 0$
    \end{definition}
    \begin{exercice}
    Discutez les stabilités A et L des trois méthodes ci-dessus. Illustrez votre réponse en codant l'équation de la chaleur et en la résolvant et montrant la solution à $t=0.1$, $N=12$, $\lambda=0.5, 1.5, 2.5$ et pour la fonction $g$ telle que définie dans l'exercice \ref{eqchal} cas 2, avec $a=\frac{1}{4}$ et $b=\frac{3}{4}$. 
    \end{exercice}
    \end{frame}
    
    %\begin{frame}
    %\includegraphics{heat_0.pdf}
    %\end{frame}
    
    \begin{frame}
    \includegraphics{heat_1.pdf}
    \end{frame}
    
    \begin{frame}
    \includegraphics{heat_xl_1.pdf}
    \end{frame}
    
    
    \begin{frame}
    \includegraphics{heat_xl_2.pdf}
    \end{frame}
    
    %\begin{frame}
    %\includegraphics{heat_3d.pdf}
    %\end{frame}
     
    
    \begin{frame}
    \begin{exercice}
      Considérez l'équation $u_t = D u_{xx} - bu$, où $D$,$b$ sont des constantes positives, $u(0,t)=u(1,t)=0$ et $u(x,0)=g(x)$.
    \begin{enumerate}
      \item Montrez que la solution exacte est $u(x,t)=\sum_{n=1}^{\infty} A_n e^{-(b+D\lambda_n^2)t} \sin(\lambda_n x)$, où $\lambda_n=n\pi$ et $A_n$ sont les coefficients de Fourier $A_n=2\int_0^1 g(x)\sin(\lambda_n x)\,\mathrm{d} x$. 
      \item Cherchez un schéma \emph{implicite} en différence finie d'erreur $\mathcal{O}(h^2)+\mathcal{O}(k)$. Exprimez le résultat sous forme matricielle. Quel est l'intérêt d'un schéma implicite dans ce contexte\ ?
      \item Déterminez la condition de stabilité.
      \item Dans le cas où $g(x)=\sin(\pi x)$, $D=0,1$ et $b=1$, imposez $\lambda=0,4$ et tracez l'erreur maximale au temps $t=1$ en fonction du nombre de pas de temps $M=4,16,64,256,1024$ (en log-log). Superposez ensuite à ce graphe les résultats obtenus pour $\lambda=0,04$.  Qu'observez-vous ? 
    \end{enumerate}
    \end{exercice}
    \end{frame}