Page MenuHomec4science
Files F62940787

neuman.tex
No OneTemporary
Actions

File Metadata

Created
Thu, May 16, 16:17

neuman.tex
View Options

%
% @file neuman.tex
%
% @brief
%
% @copyright
% Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
% SPC (Swiss Plasma Center)
%
% spclibs is free software: you can redistribute it and/or modify it under
% the terms of the GNU Lesser General Public License as published by the Free
% Software Foundation, either version 3 of the License, or (at your option)
% any later version.
%
% spclibs is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with this program. If not, see <https://www.gnu.org/licenses/>.
%
% @authors
% (in alphabetical order)
% @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
%
\documentclass[a4paper]{article}
\usepackage{amsmath}
\title{\tt Some Notes on Boundary Conditions}
\author{Trach-Minh Tran}
\date{March 2012}
\begin{document}
\maketitle
\section{Neumann BC as an {essential} BC}
The original equation:
\begin{equation}
\mathbf{A \cdot u} = \mathbf{b}
\end{equation}
with the Neumann BC (1D case):
\begin{equation}
\alpha u_1 + \beta u_2 = c.
\end{equation}
From Eq.(20) of \cite{BSPLINES}:
\begin{equation}
\beta = -\alpha =\frac{p}{\Delta_1}
\end{equation}
where $p$ is the degree of spline and $\Delta_1$ is the lenght of the
first insterval.
Transformation $(u_1, \ldots, u_n) \Rightarrow (\hat u_1, \ldots, \hat
u_n)$ defined by
\begin{equation}
\begin{array}{ccc}
\alpha u_1 + \beta u_2 = \hat u_1 & & u_1 = \frac{1}{\alpha}\hat u_1 - \frac{\beta}{\alpha}\hat u_2 \\
u_2 = \hat u_2 & & u_2 = \hat u_2 \\
\vdots & \Longrightarrow & \vdots \\
u_N = \hat u_N & & u_N = \hat u_N.
\end{array}
\end{equation}
The original Neumann BC becomes now a \emph{inhomogeneous Dirichlet}
BC on $\mathbf{\hat u}$:
\begin{equation}
\hat u_1 = c.
\end{equation}
The transformed linear system can be written as:
\begin{equation}
\mathbf{(U^T\cdot A \cdot U)\cdot \hat u} = \mathbf{U^T\cdot b},
\end{equation}
where $\mathbf{U}$ is given by
\begin{equation}
\mathbf{U} =
\left(\begin{matrix}
\frac{1}{\alpha} & -\frac{\beta}{\alpha} & \dots & 0 \\
0 & 1 & \dots & 0 \\
& & \ddots& \vdots \\
0 & 0 & \dots & 1
\end{matrix}\right)
\end{equation}
Thus, all the symmetry, hermiticity or positivity properties of the
original matrix are preserved with this matrix transformation!
\section{Neumann BC as a \emph{natural} BC}
Multiplying the 1D Sturm-Liouville equation (see section 1.1.1 of
\cite{SOLVERS}) by spline $\Lambda_j(x)$ and integrating by parts, we
obtain the following boundary terms:
\begin{equation}
-\Lambda_j(L) C_1(L) \phi'(L) + \Lambda_j(0) C_1(0) \phi'(0)
\end{equation}
To impose $\phi'(0) = a$ and noting that $\Lambda_j(0)=\delta_{j1}$,
you only need to add $[-aC_1(0)]$ to the first element of the
RHS. Likewise, for the BC $\phi'(L) = b$ you only need to add
$[bC_1(L)]$ to the last element of the RHS. No matrix manipulation (as
for the \emph{essential} BC) is required! Notice that if $a$ or $b$ is
zero, nothing needs to be done to impose these BC. That's the reason why
it is called \emph{natural} BC!
A subtle point to be noted here is that using \emph{natural} BC,
$\phi'(0)$ \emph{is not} exaclty equal to $a$, althought it should
converge to $a$ as $(\Delta x)^p$ where $p$ is the spline degree,
while using the \emph{essential} BC, $\phi'(0)=a$ is \emph{exact}!
\section{Diffusion Equation using second order time implicit method}
Let rewrite Eq.(74) of your notes in vector form and replace the
unkowns $n$ by $f$:
\begin{equation}
\mathbf{B} \frac{d \mathbf{f}}{dt} = \mathbf{M\cdot f}.
\end{equation}
Using a \emph{second order time centered} discretization,
\begin{equation}
\begin{split}
\mathbf{B} \left(\frac{\mathbf{f}^{n+1}-\mathbf{f}^{n}}{\Delta t}\right) &=
\mathbf{M} \left(\frac{\mathbf{f}^{n+1}+\mathbf{f}^{n}}{2}\right) \\
\Rightarrow &
\left(\mathbf{B} -\frac{\Delta t}{2}
\mathbf{M}\right)\mathbf{f}^{n+1} =
\left(\mathbf{B} +\frac{\Delta t}{2}
\mathbf{M}\right)\mathbf{f}^{n}
\end{split}
\end{equation}
\emph{Essential} BC has to be imposed on the matrix
\begin{equation}
\mathbf{B} -\frac{\Delta t}{2} \mathbf{M}
\end{equation}
while \emph{natural} BC is introduced while deriving the weak form
leading to the matrix $M$.
This method is \emph{unconditionnaly stable} and second order in
time. When linear splines are used for the space discretization, this
scheme is similar to the well-known \emph{Cranck-Nicolson} (see for
example Wikipedia) discretization for parabolic PDE.
\begin{thebibliography}{99}
\bibitem{BSPLINES} {\tt BSPLINES} Reference Guide.
\bibitem{SOLVERS} {\tt The SOLVERS in BSPLINES} Reference Guide.
\end{thebibliography}
\end{document}

Event Timeline

c4science · Help