readme.tex

解决麦克斯韦的matlab源文件

\documentclass[12pt]{article}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{amscd}
\usepackage{latexsym}
\usepackage{citesort}
\usepackage[dvips]{graphicx} 
\usepackage{psfrag}
\usepackage{a4}
\usepackage{sober}
%\usepackage[draft]{showkeys}

\setlength{\topmargin}{-1.5cm} 
\setlength{\oddsidemargin}{0.1cm}
\setlength{\textwidth}{16.0cm} 
\setlength{\textheight}{24.0cm}
\setlength{\parskip}{6pt} 
\setlength{\parindent}{0pt} 
\frenchspacing
\sloppy

\input{symbols}
\input{symb_ops}
\input{symb_spaces}
\input{symb_fem}

\begin{document}
\DeclareGraphicsRule{.eps.bz2}{eps}{.eps.bb}{`bunzip2 -c #1}
\DeclareGraphicsRule{.eps.bz}{eps}{.eps.bb}{`bunzip2 -c #1}
\DeclareGraphicsRule{.eps.gz}{eps}{.eps.bb}{`gunzip -c #1}
(C) D. Arnold, R.Hiptmair, \today \hfill Not for dissemination
\vspace{1.5cm}

\centerline{\Large Discretizing transient Maxwell's equations}
\vspace{1cm}

\section{Continuous problem}

On a bounded domain $\Omega\subset\bbR^{2}$ we consider the 
electric wave equation for the electric field $\VE=\VE(\Vx,t)$,
$\Vx\in\Omega$, $t\in[0;T]$,
\begin{gather}
  \label{eq:1}
  \tfrac{d^{2}}{dt^{2}}\VE + \plcurl\scurl\VE = 0\;,\\
  \VE(.,t)\times\Vn = 0 \quad\forall t\;,\\
  \VE(\Vx,0) = \VE_{0}\in\nHcurl\quad,\quad \tfrac{d}{dt}\VE(\Vx,0) = 0\;,
\end{gather}
with the differential operators
\begin{gather*}
  \plcurl u = 
  \begin{pmatrix}
    -\frac{\partial u}{\partial y} & \frac{\partial u}{\partial x}
  \end{pmatrix}^{T}\quad,\quad
  \scurl\Vu = \frac{\partial u^{1}}{\partial y} - \frac{\partial u^{2}}{\partial x}\;.
\end{gather*}
It can be derived from the $z$-invariant Maxwell's equations assuming homogeneous,
isotropic materials inside $\Omega$. In the sequel we assume that $\Omega$ is a
polygon. Moreoever, we take for granted that
\begin{gather}
  \label{df}
  \Div\VE_{0} = 0\;.
\end{gather}

The variational formulation of \eqref{eq:1} reads: seek $\VE\in C^{2}([0;T],\nHcurl)$
such that
\begin{gather}
  \label{eq:2}
  \SPLzweiv{\tfrac{d^{2}}{dt^{2}}\VE}{\Vv} + \SPLzweiv{\scurl\VE}{\scurl\Vv} = 0
  \quad\forall \Vv\in\nHcurl\;,\\
  \VE(\Vx,0) = \VE_{0}\quad,\quad \tfrac{d}{dt}\VE(\Vx,0) = 0\;.
\end{gather}

\section{Spatial discretization}

We opt for a conforming spatial finite element discretization of \eqref{eq:2}.  To
that end we equip $\Omega$ with a triangular mesh $\Ct_{h}=\{T_{i}\}_{i}$.  The
following finite element spaces will be relevant
\begin{align*}
  \Cs_{h} & 
  := \{v\in C^{0}(\Omega)\cap\nHeins,\; v_{|T}\in\Cp_{1}(T)\,\forall T\in\Ct_{h}\}\;,\\
  \Cev_{h} & := \text{lowest order Whitney 1-forms } \subset \nHcurl\;,\\
  \Cnv_{h} & := \{\Vv\in (C^{0}(\Omega))^{2}\cap\nHcurl,\; \Vv_{|T}\in(\Cp_{1}(T))^{2}
  \,\forall T\in\Ct_{h}\}\;.
\end{align*}

\section{Regularization}

In principle, condition \eqref{df} ensures that $\Div\VE(.,t)=0$ for all times.
However, a slight perturbation of the initial value might lead to growing $\scurl$-free
components in $\VE(.,t)$ that may eventually swamp the physically meaningful solution.
A remedy is offered by regularization. 

\subsection{Continuous regularization}

This strategy swaps \eqref{eq:2} for the completely equivalent variational
problem: seek $\VE\in C^{2}([0;T],\nHcurl\cap \Hdiv)$ such that for all
$\Vv\in\nHcurl\cap \Hdiv$
\begin{gather}
  \label{reg}
  \SPLzweiv{\tfrac{d^{2}}{dt^{2}}\VE}{\Vv} + \SPLzweiv{\scurl\VE}{\scurl\Vv}
  + \SPLzwei{\Div\VE}{\Div\Vv} = 0 \;,\\
   \VE(\Vx,0) = \VE_{0}\quad,\quad \tfrac{d}{dt}\VE(\Vx,0) = 0\;.
\end{gather}
This variational problem can be discretized in $\Cnv_{h}$ provided that we
use a projection of $\VE_{0}$ as initial guess. We arrive at the linear
system of equations
\begin{gather}
  \label{asdfger}
  \hat{\Msf}\tfrac{d^{2}}{dt^{2}}\vec{\VE} + \hat{\Csf}\vec{\VE}
  + \hat{\Rsf}\vec{\VE} = 0\;.
\end{gather}
Here $\hat{\Msf}\in\bbR^{N,N}$, $N:=\Dim =Cnv_{h}$, is the mass matrix
with respect to $\Cnv_{h}$, $\hat{\Csf}$ the discrete counterpart of
$\plcurl\scurl$, and $\hat{\Rsf}$ corresponds to the discrete regularization
term. 

However, \eqref{asdfger} might not be a valid discretization of \eqref{reg}, because
$\Heinsv$ fails to be dense in $\nHcurl\cap \Hdiv$ whenever $\Omega$ features
reentrant corners. It is our objective to demonstrate this effect in a 
numerical experiment. 

\subsection{Discrete regularization}

The variational problem \eqref{reg} is not amenable to a discretization by means
of edge elements. The alternative discrete regularization procedure is based  
on the observation that, thanks to $\grad\Cs_{h}\subset\Cev_{h}$, we have
\begin{gather*}
  \SPLzweiv{\VE_{h}(t)}{\grad v_{h}} = 0\quad \forall v_{h} \in \Cs_{h}\;,
\end{gather*}
if $\VE_{h}\in\Cev_{h}$ is the solution of \eqref{eq:2} after Galerkin discretization
in space using $\Cev_{h}$. Therefore, $\VE_{h}(t)$ can also be obtained as the solution
of the saddle point problem: seek $\VE_{h} \in C^{2}([0;T],\Cev_{h})$, $p_{h}\in \Cs_{h}$,
such that
\begin{gather}
  \label{dreg}
  \begin{array}[c]{ccccccll}
    \SPLzweiv{\frac{d^{2}}{dt^{2}}\VE_{h}}{\Vv_{h}}&+&
    \SPLzweiv{\scurl\VE_{h}}{\scurl\Vv_{h}}&+&\SPLzweiv{\Vv_{h}}{\grad p_{h}}&=& 0\,
    &\forall\Vv_{h}\in\Cev_{h}\;,\\
    &&\SPLzweiv{\VE_{h}}{\grad q_{h}}&-&d(p_{h},q_{h})&=&0\,&\forall q_{h}\in\Cs_{h}\;,
  \end{array}
\end{gather}
where $d(\cdot,\cdot)$ is an arbitrary positive definite bilinear form on $\Cs_{h}$.
To see this, note that $p_{h}=0$. We have decided to use the lumped $\Lzwei$ inner product
on $\Cs_{h}$ in order to introduce $d(\cdot,\cdot)$. Then \eqref{dreg} gives rise
to the following linear system of equations
\begin{gather}
  \label{mayreg}
   \begin{array}[c]{ccccccl}
     \Msf \frac{d^{2}}{dt^{2}}\vec{\VE} &+& \Csf\vec{\VE}&+& \Gsf\vec{p}&=&0\;,\\
     && \Gsf^{T}\vec{\VE}&-& \Dsf\vec{p}&=& 0\;,
   \end{array}
 \end{gather}
where $\Dsf$ is diagonal, $\Msf$ the edge element mass matrix, $\Csf$ the edge
element $\plcurl\scurl$-matrix, and $\Gsf$ the discrete gradient. The linear
system \eqref{mayreg} can be recast as
\begin{gather}
  \label{sdfg}
  \Msf \frac{d^{2}}{dt^{2}}\vec{\VE} + (\Csf+\Gsf\Dsf^{-1}\Gsf^{T})\vec{\VE} = 0\;.
\end{gather}
It goes without saying that suitable projected initial values have to be 
specified. 

\section{Timestepping}
\label{sec:timestepping}

The second derivative in time will be discretized by means of the explicit trapezoidal rule
(leap-frog), which is a two-step method
\begin{gather*}
  \frac{d^{2}}{dt^{2}}\VE_{h} \approx \dfrac{\vec{\VE}^{n+1}-2\vec{\VE}^{n}+\vec{\VE}^{n-1}}{\tau^{2}}\;,
\end{gather*}
where $\tau$ is the size of the timestep. It has to be sufficiently small to satisfy
a CFL-condition. In ech timestep we determine $\vec{\VE}^{n+1}$ from
$\vec{\VE}^{n}$ and $\vec{\VE}^{n-1}$ according to
\begin{align*}
  \text{For }\Cev_{h}:\quad & \vec{\VE}^{n+1} = 2\vec{\VE}^{n}-\vec{\VE}^{n-1} + \tau^{2}
  \Msf^{-1}(\Csf+\Gsf\Dsf^{-1}\Gsf^{T})\vec{\VE}^{n}\;,\\
  \text{For }\Cnv_{h}:\quad &
  \vec{\VE}^{n+1} = 2\vec{\VE}^{n}-\vec{\VE}^{n-1} + \tau^{2}
  \hat{\Msf}^{-1}(\hat{\Csf}+\hat{\Rsf})\vec{\VE}^{n}\;.
\end{align*}
In addition, $\vec{\VE}^{0}$ are the nodal values of a suitable projection of 
$\VE_{0}$ onto $\Cev_{h}$ and $\Cnv_{h}$, respectively, and $\vec{\VE}^{-1}=\vec{\VE}^{1}$.
Therefore, the first step boils down to
\begin{gather*}
  \vec{\VE}^{1} = \vec{\VE}^{0} +
  \tfrac{1}{2}\tau^{2}\Msf^{-1}(\Csf+\Gsf\Dsf^{-1}\Gsf^{T})\vec{\VE}^{0}\;.
\end{gather*}
A crucial issue is the computation of the discrete initial guess $\vec{\VE}^{0}$.  It
is supposed that a divergence-free $\VE_{0}$ is given as an analytic expression.
\begin{itemize}
\item In the case of vertex based discretization in $\Cnv_{h}$ we take
  $\vec{\VE}^{0} := I_{h}\VE$, where $I_{h}$ is the interpolation in
  the vertices of the mesh.
\item In the case of edge elements it may be important that $\vec{\VE}^{0}$ is
  discretely divergence free, that is
  \begin{gather*}
    \SPLzweiv{\VE_{h}^{0}}{\grad\phi_{h}} = 0 \quad\forall \phi_{h}\in \Cs_{1}\;.
  \end{gather*}
  An option is to determine a solution of the saddle point problem: seek
  $\VE_{h}^{0}\in\Cnv_{h}$, $\Vu_{h}\in\Cnv$ such that for all $\Vv_{h}\in\Cnv_{h}$
  and $\Vw_{h}\Cnv_{h}$
  \begin{gather}
    \label{eq:3}
    \begin{array}[c]{ccccl}
      \SPLzweiv{\VE_{h}^{0}}{\Vv_{h}}&+&\SPLzwei{\scurl\Vu_{h}}{\scurl\Vv_{h}}&=&
      \SPLzweiv{\VE_{0}}{\Vv_{h}}\;,\\
      \SPLzweiv{\scurl\VE_{h}^{0}}{\scurl\Vw_{h}}&&&=&
      \SPLzweiv{\scurl\VE^{0}}{\scurl\Vw_{h}}\;.
    \end{array}
  \end{gather}
  Let $\Pi_{h}$ stand for the local edge element interpolation operator. In two
  dimensions the commuting diagram property can be stated as
  \begin{gather*}
    \scurl\circ\Pi_{h} = Q_{h}\circ \scurl\;,
  \end{gather*}
  where $Q_{h}$ is the $\Lzwei$-orthogonal projection onto the space of
  $\Ct_{h}$-piecewise constant functions. Thus, $\scurl(\VE_{h}^{0}-\Pi_{h}\VE_{0})=0$
  holds for the solution $\VE_{h}^{0}$ of \eqref{eq:3}. Assuming simple topology
  of $\Omega$ this implies
  \begin{gather*}
    \exists\phi_{h}\in\Cs_{1}\;:\quad\VE_{h}^{0}-\Pi_{h}\VE_{0} = \grad\phi_{h}\;.
  \end{gather*}
  It satisfies, for all $\psi_{h}\in\Cs_{1}$,
  \begin{gather*}
    \SPLzweiv{\grad\phi_{h}}{\grad\psi_{h}} =
    \SPLzweiv{\VE_{0}-\Pi_{h}\VE_{0}}{\grad\psi_{h}} = 
    -\SPLzweiv{\Pi_{h}\VE_{0}}{\grad\psi_{h}}\;.
  \end{gather*}
  This amounts to the linear system of equations
  \begin{gather*}
    \Gsf^{T}\Msf\Gsf\vec{\phi} = \Gsf^{T}\Msf\,\vec{\Pi_{h}\VE_{0}}\;,
  \end{gather*}
  so that we end up with
  \begin{gather}
    \label{getinti}
    \vec{\VE}^{0} = 
    (Id - \Gsf(\Gsf^{T}\Msf\Gsf)^{-1}\Gsf^{T}\Msf)\vec{\Pi_{h}\VE_{0}}\;.
  \end{gather}
\end{itemize}

\newcommand{\pfrac}[2]{\frac{\partial{#1}}{\partial{#2}}}

For the L-shaped domain $\Omega=]-1,1[\times]0;1[\cup[0;1]\times]-1;0[$ we 
choose an initial field that contains a singular contribution at the 
reentrant corner at $(0,0)$. We start with a singular function for
the Laplacian (in polar coordinates), $\omega=\frac{3}{2}\pi$
\begin{gather*}
  u(r,\varphi) := r^{\pi/\omega}\sin(\frac{\pi}{\omega}\,\varphi)\quad
  r\geq 0, 0<\varphi<\frac{3}{2}\pi\;.
\end{gather*}
The associated singular vector field reads
\begin{gather*}
  \grad u(r,\varphi) = \pfrac{u}{r}\cdot\vec{e}_{r} + \frac{1}{r}\pfrac{u}{\varphi}
  \cdot\vec{e}_{\varphi}\;.
\end{gather*}
Then, we introduce 
\begin{gather*}
  p(r,\varphi) := r^{\pi/\omega}\cos(\frac{\pi}{\omega}\,\varphi)\;,
\end{gather*}
and see
\begin{gather*}
  \left.
  \begin{array}[c]{rcl}
      \pfrac{p}{r} &=& \frac{1}{r}\pfrac{u}{\varphi}\\
      -\frac{1}{r}\pfrac{p}{\varphi} &=& \pfrac{u}{r}
  \end{array}\right\}
\quad
\Rightarrow
\quad
\plcurl p := 
\pfrac{p}{r}\cdot\vec{e}_{\varphi} - \frac{1}{r}\pfrac{p}{\varphi}\cdot\vec{e}_{r}
= \grad u\;.
\end{gather*}
Next, we have to resort to a cut-off function: pick
\begin{gather*}
  g(t) :=
  \begin{cases}
    1 & \text{ if } 0\leq |t| \leq \tfrac{1}{2}\;,\\
    \sin^{2}(\pi\,t) & \text{ if } \tfrac{1}{2}\leq |t| \leq 1
  \end{cases}\quad\Rightarrow\quad
  g'(t) = 
  \begin{cases}
    0 & \text{ if } 0\leq |t| \leq \tfrac{1}{2}\;,\\
    \pi\sin(2\pi\,t) & \text{ if } \tfrac{1}{2}\leq |t| \leq 1
  \end{cases}\;.
\end{gather*}
We write $f(x,y) = g(x)g(y)$ and define
\begin{gather}
  \label{En}
  \VE^{0}(x,y) := \plcurl (p(r,\varphi)\cdot f(x,y)) = 
  f\,\plcurl p + p\,\plcurl f\;.
\end{gather}
Note that
\begin{gather*}
  \plcurl f(x,y) = 
  \begin{pmatrix}
    -g(x)g'(y) \\ g'(x)g(y)
  \end{pmatrix} = 
  \begin{cases}
    \begin{pmatrix}
      0 \\ 0
    \end{pmatrix}
    & \text{ if } 0\leq|x|,|y|\leq\tfrac{1}{2}\;,\\
    \begin{pmatrix}
      0 \\
      \pi\sin(2\pi\,x)
    \end{pmatrix} &
    \text{ if } 0\leq|y|\leq\tfrac{1}{2},|x|\geq\tfrac{1}{2}\;,\\
    \begin{pmatrix}
      -\pi\sin(2\pi\,y) \\
      0 
    \end{pmatrix} & 
    \text{ if } 0\leq|x|\leq\tfrac{1}{2},|y|\geq\tfrac{1}{2}\;,\\
    \begin{pmatrix}
      -\pi\sin^{2}(\pi x)\sin(2\pi y) \\
      \pi\sin(2\pi x)\sin^{2}(\pi y)
    \end{pmatrix} &
    \text{ if } \tfrac{1}{2}\leq |x|,|y| \leq 1
  \end{cases}\;.
\end{gather*}

\section{Element matrices}

Let $T$ denote a single triangle with vertices $\Va_i=(x_i,y_i)$, $i=1,2,3$. Let
$\lambda_i$ denote the barycentric coordinate functions, i.e., the linear functions
with $\lambda_i(\Va_j)=\delta_{ij}$. If the vertices are arranged in counterclockwise
fashion, they have the representation
\begin{gather*}
  \lambda_{1}(\Vx) = \frac{1}{2|T|}(\Vx-
  \begin{pmatrix}
    x_{2} \\ y_{2}
  \end{pmatrix}
  )\cdot
  \begin{pmatrix}
    y_{2}-y_{3}\\
    x_{3}-x_{2}
  \end{pmatrix}\;,\\
  \lambda_{2}(\Vx) = \frac{1}{2|T|}(\Vx-
  \begin{pmatrix}
    x_{3} \\ y_{3}
  \end{pmatrix}
  )\cdot
  \begin{pmatrix}
    y_{3}-y_{1}\\
    x_{1}-x_{3}
  \end{pmatrix}\;,\\
  \lambda_{3}(\Vx) = \frac{1}{2|T|}(\Vx-
  \begin{pmatrix}
    x_{1} \\ y_{1}
  \end{pmatrix}
  )\cdot
  \begin{pmatrix}
    y_{1}-y_{2}\\
    x_{2}-x_{1}
  \end{pmatrix}\;.
\end{gather*}
Their gradients are constant and read
\begin{gather*}
  \grad\lambda_{1} = \frac{1}{2|T|}\begin{pmatrix}
    y_{2}-y_{3}\\
    x_{3}-x_{2}
  \end{pmatrix}\quad,\quad
  \grad\lambda_{2} = 
  \frac{1}{2|T|}\begin{pmatrix}
    y_{3}-y_{1}\\
    x_{1}-x_{3}
  \end{pmatrix}\quad,\quad
  \grad\lambda_{3} = 
  \frac{1}{2|T|}\begin{pmatrix}
    y_{1}-y_{2}\\
    x_{2}-x_{1}
  \end{pmatrix}\;.
\end{gather*}

\begin{figure}[!htbp]
  \centering
  \psfrag{a1}{$\Va_{1} = \binom{x_{1}}{y_{1}}$}
  \psfrag{a2}{$\Va_{2} = \binom{x_{2}}{y_{2}}$}
  \psfrag{a3}{$\Va_{3} = \binom{x_{3}}{y_{3}}$}
  \includegraphics[width=0.6\textwidth]{tri.eps}
  \caption{Generic triangular element}
  \label{fig:tri}
\end{figure}

Given a bilinear form $b_{T}$ and
local basis functions $\phi_{1},\ldots,\phi_{m}$ we aim to compute the element
matrix $\Ssf_{T}\in\bbR^{m,m}$ with
\begin{equation*}
\Ssf_{ij}=b_T(\phi_j,\phi_i).
\end{equation*}

\subsection{Space $\Cs_{h}$}

In this space we only need a lumped mass matrix, obtained from the local
bilinear form
\begin{gather*}
  b_{T}(\phi,\psi) = \frac{1}{3}|T|\sum\limits_{j=1}^{3} \phi(\Va_{j})\psi(\Va_{j})\;.
\end{gather*}
The lumped local mass matrix is diagonal and reads $\Dsf_{T} = \frac{1}{3}|T|\, I_{3}$.
 
\subsection{Space $\Cnv_{h}$}

The basis functions are (in the order we shall take) the 2-vector-valued functions
\begin{equation*}
(\lambda_1,0),\quad(0,\lambda_1),\quad(\lambda_2,0),
\quad(0,\lambda_2),\quad(\lambda_3,0),\quad(0,\lambda_3)\;.
\end{equation*}

\subsubsection{Local mass matrix}

The local bilinear form agrees with the $\xLzwei{T}$ inner product.
From the scalar case, we see that the local mass matrix is
\begin{equation*}