readme.tex

解决麦克斯韦的matlab源文件

\hat{\Msf}_{T}=\frac{|T|}{12}\begin{pmatrix} 2&0&1&0&1&0\\0&2&0&1&0&1\\
 1&0&2&0&1&0\\0&1&0&2&0&1\\1&0&1&0&2&0\\0&1&0&1&0&2\end{pmatrix}.
\end{equation*}

\subsection{Local $\plcurl\scurl$ matrix}

\newcommand*{\pd}[2]{\frac{\partial #1}{\partial #2}}
The local bilinear form is
\begin{equation*}
b_T(\phibf,\psibf)=\int_T\scurl\psibf\,\scurl\phibf\,d\Vx\;,
\end{equation*}
where $\scurl\phibf=\frac{\partial\phibf^1}{\partial y}-\frac{\partial\phibf^2}{\partial x}$.
We get $\hat{\Csf}_{T}=|T|\Ksf^{T}\Ksf$, where $\Ksf$ is the $1\times6$ matrix whose
entries are $\scurl\phibf_i$, i.e.,
\begin{align*}
\Ksf&=\begin{pmatrix}\displaystyle\pd{\lambda_1}y,-\pd{\lambda_1}x,\pd{\lambda_2}y,-\pd{\lambda_2}x,\pd{\lambda_3}y,-\pd{\lambda_3}x\end{pmatrix}
\\
&=\frac{1}{2|T|}
\begin{pmatrix}
x_3-x_2,y_3-y_2,x_1-x_3,y_1-y_3,x_2-x_1,y_2-y_1
\end{pmatrix}.
\end{align*}

\subsection{Local $-\grad\Div$ matrix}

The local bilinear form ist
\begin{gather*}
  b_{T}(\phibf,\psibf)= \int_{T}\Div\phibf\,\Div\psibf\,d \Vx
\end{gather*}
We form the row vector
\begin{align*}
  \Lsf & := 
  \begin{pmatrix}
    \Div\phibf_{1} & 
    \Div\phibf_{2} & 
    \Div\phibf_{3} & 
    \Div\phibf_{4} & 
    \Div\phibf_{5} & 
    \Div\phibf_{6}
  \end{pmatrix} \\
  & = 
  \frac{1}{2|T|}
  \begin{pmatrix}
    (y_{2}-y_{3}) & 
    (x_{3}-x_{2}) & 
    (y_{3}-y_{1}) & 
    (x_{1}-x_{3}) & 
    (y_{1}-y_{2}) & 
    (x_{2}-x_{1})
  \end{pmatrix}\;.
\end{align*}
It permits us to express $\hat{\Rsf}_{T} = |T|\Lsf^{T}\Lsf$.

\subsection{Space $\Cev_{h}$}

We denote by $e_{i}$ the edge of $T$ opposite to the vertex $\Va_{i}$.  It is
supposed to be directed from $\Va_{i+1}$ to $\Va_{i+2}$, where all indices have to be
read $\mod 3 + 1$. The local edge element basis functions are given by
\begin{align*}
  \phibf_{1}(\Vx) & := 
  \lambda_{2}(\Vx)\grad\lambda_{3} - \lambda_{3}(\Vx)\grad\lambda_{2}\;,\\
  \phibf_{2}(\Vx) & := 
  \lambda_{3}(\Vx)\grad\lambda_{1} - \lambda_{1}(\Vx)\grad\lambda_{3}\;,\\
  \phibf_{3}(\Vx) & := 
  \lambda_{1}(\Vx)\grad\lambda_{2} - \lambda_{2}(\Vx)\grad\lambda_{1}\;.
\end{align*}
These local vector fields satisfy
\begin{gather*}
  \int\limits_{e_{i}}\phibf_{j}\cdot d\vec{s} = \delta_{ij}\;.
\end{gather*}
The $\scurl$ of these basis functions is constant on $T$.
By means of Gau\ss' theorem, we can compute 
\begin{gather*}
  \scurl\phibf_{i} = \frac{1}{|T|}\int\limits_{T}\scurl\phibf_{i}\,d\Vx
  = -\frac{1}{|T|} \int\limits_{\partial T}\phibf\cdot \,d\vec{s} = 
  - \frac{1}{|T|}\;.
\end{gather*}
At the barycenter $\Vc_{T}$ of $T$ the basis functions evaluate to
\begin{align*}
  \phibf_{1}(\Vc_{T}) & = \frac{1}{3}(\grad\lambda_{3}-\grad\lambda_{2})\;,\\
  \phibf_{2}(\Vc_{T}) & = \frac{1}{3}(\grad\lambda_{1}-\grad\lambda_{3})\;,\\
  \phibf_{3}(\Vc_{T}) & = \frac{1}{3}(\grad\lambda_{2}-\grad\lambda_{1})\;.
\end{align*}

\subsubsection{Mass matrix}

From the above representations of the local basis functions we easily obtain
their values at the midpoint $\Vm_{i}$ of edge $e_{i}$:
\begin{gather*}
  \begin{array}[c]{c}
    \phibf_{1}(\Vm_{1}) = \frac{1}{4|T|}
    \begin{pmatrix}
      2y_{1} - y_{3} - y_{2} \\
      -2x_{1} + x_{3} + x_{2}
    \end{pmatrix} ,\,
    \phibf_{1}(\Vm_{2}) = \frac{1}{4|T|}
    \begin{pmatrix}
      y_{1} - y_{3} \\
      x_{3} - x_{1}
    \end{pmatrix} ,\,
    \phibf_{1}(\Vm_{3}) = \frac{1}{4|T|}
    \begin{pmatrix}
      y_{1} - y_{2} \\
      x_{2} - x_{1}
    \end{pmatrix}\;, \\
    \phibf_{2}(\Vm_{1}) = \frac{1}{4|T|}
    \begin{pmatrix}
      y_{2}-y_{3} \\
      x_{3}-x_{2}
    \end{pmatrix} ,\, 
    \phibf_{2}(\Vm_{2}) = \frac{1}{4|T|}
    \begin{pmatrix}
      2y_{2}-y_{3}-y_{1}\\
      -2x_{2}+x_{3}+x_{1}
    \end{pmatrix} ,\, 
    \phibf_{2}(\Vm_{3}) = \frac{1}{4|T|}
    \begin{pmatrix}
      y_{2}-y_{1} \\
      x_{1}-x_{2}
    \end{pmatrix}\;, \\
    \phibf_{3}(\Vm_{1}) = \frac{1}{4|T|}
    \begin{pmatrix}
      y_{3}-y_{2}\\
      x_{2}-x_{3}
    \end{pmatrix} ,\,
    \phibf_{3}(\Vm_{2}) = \frac{1}{4|T|}
    \begin{pmatrix}
      y_{3}-y_{1} \\
      x_{1}-x_{3}
    \end{pmatrix} ,\,
    \phibf_{3}(\Vm_{3}) = \frac{1}{4|T|}
    \begin{pmatrix}
      2y_{3}-y_{1}-y_{2}\\
      -2x_{3}+x_{1}+x_{2}
    \end{pmatrix}\;.
  \end{array}
\end{gather*}
The local mass matrix is associated with the bilinear form
\begin{gather*}
  b_{T}(\phibf,\psibf) := \int\nolimits_{T}\phibf\cdot\psibf\,d\Vx\;.
\end{gather*}
As the integrand is a quadratic polynomial, the entries of the local mass matrix can
be evaluated exactly using a midpoint based quadrature formula
\begin{gather*}
  b_{T}(\phibf,\psibf) \approx
  \frac{1}{3}|T|\sum\limits_{j=1}^{3}\phibf(\Vm_{j})\cdot\psibf(\Vm_{j})\;.
\end{gather*}
In order to obtain a concise expression for the local mass matrix we introduce
\begin{gather*}
  \Psf := 
  \begin{pmatrix}
    2y_{1} - y_{3} - y_{2}  & y_{2}-y_{3}         & y_{3}-y_{2}   \\
    -2x_{1} + x_{3} + x_{2} & x_{3}-x_{2}         & x_{2}-x_{3}   \\
    y_{1} - y_{3}           & 2y_{2}-y_{3}-y_{1}  & y_{3}-y_{1}   \\
    x_{3} - x_{1}           & -2x_{2}+x_{3}+x_{1} & x_{1}-x_{3}   \\
    y_{1} - y_{2}           & y_{2}-y_{1}         & 2y_{3}-y_{1}-y_{2} \\
    x_{2} - x_{1}           & x_{1}-x_{2}         & -2x_{3}+x_{1}+x_{2} \\
  \end{pmatrix}\;.
\end{gather*}
Then the $3\times3$ local mass matrix for lowest order edge elements
can be written as
\begin{gather*}
  \Msf_{T} = \frac{1}{42|T|}\Psf^{T}\Psf\;.
\end{gather*}

\subsubsection{$\plcurl\scurl$-Matrix}

The relevant local bilinear form is
\begin{gather*}
  b_{T}(\phibf,\psibf) = \int\nolimits_{T}\scurl\phibf\,\scurl\psibf\,d\Vx\;.
\end{gather*}
Due to the constant $\scurl$s of the local basis functions, we can simply
write $\Ssf_{T} := |T|\Qsf^{T}\Qsf$, where 
\begin{gather*}
  \Csf := \frac{1}{|T|}
  \begin{pmatrix}
    1 & 1 & 1
  \end{pmatrix}
  \in\bbR^{1,3}\;.
\end{gather*}

\subsubsection{Discrete gradient}

From topological considerations the discrete gradient matrix $\Gsf$ emerges
as the incidence matrix between vertices and oriented edges: if $e=[\Vv_{1},\Vv_{2}]$ is an edge of the global mesh, then 
\begin{gather*}
  \Gsf_{e\Vv_{1}} = -1\quad,\quad \Gsf_{e\Vv_{2}} = 1\;.
\end{gather*}

\section{MATLAB implementation}

\begin{itemize}
\item \texttt{leapfrog.m}: Main entry point and timestepping. 
\item \texttt{getinitval.m}: Computation of discrete initial values for both
  edge elements and nodal elements.
\item \texttt{Emats.m}: Routine providing mass matrix and stiffness matrix
  for edge elements
\item \texttt{Nmats.m}: Procedure for computing mass matrix $\hat{\Msf}$/stiffness
  matrix $\hat{A}$ for finite element space $\Cnv_{h}$
\item \texttt{run.m}: Main MATLAB script running the numerical experiments
\texttt{render.m}: Graphical display of dat recorded during timestepping
\end{itemize}

For the other M-files we refer to the documentation included into
the source code.

\section{Numerical experiments}

The numerical experiments were carried out on two domain, the unit square
$\Omega=]0;1[^{2}$, and the L-shaped domain $\Omega:=]-1,1[\times]0;1[\cup
]0;1[\times]-1;0[$. Both domains were equipped with triangular meshes, which
are truly unstructured for the L-shaped domain, see figure \ref{fig:lmesh}.
On the square we chose the regular triangulation created by dividing each
square of a tensor-product grid into two triangles.

\begin{figure}[!htbp]
  \centering
  \begin{minipage}[c]{0.33\textwidth}
    \begin{center}
      \includegraphics[width=\textwidth,totalheight=1.1\textwidth]{Lshape0.eps.bz2}
    \end{center}
  \end{minipage}%%
  \begin{minipage}[c]{0.33\textwidth}
    \begin{center}
      \includegraphics[width=\textwidth,totalheight=1.1\textwidth]{Lshape1.eps.bz2}
    \end{center}
  \end{minipage}%%
  \begin{minipage}[c]{0.33\textwidth}
    \begin{center}
      \includegraphics[width=\textwidth,totalheight=1.1\textwidth]{Lshape2.eps.bz2}
    \end{center}
  \end{minipage}%%
  \caption{Unstructured meshes on L-shaped domain}
  \label{fig:lmesh}
\end{figure}

The initial values $\vec{\VE}^{0}$ were chosen as explained in
Sect.~\ref{sec:timestepping} with $\omega=\frac{3}{2}\pi$ for the L-shaped domain,
and $\omega=\frac{1}{2}\pi$ for the square. The timestep was chosen
as $\tau=0.01$, $\tau=0.0025$, $\tau=0.0006125$ on the coarse, intermediate,
and fine meshes, respectively.

\subsection{Discrete evolutions with regularization}

The regularized formulation (\ref{reg}) was used in conjunction with 
globally continuous vectorial finite elements. For edge elements the
experiment were based on the discretely regularized formulation
(\ref{sdfg}). 

\begin{figure}[htbp]
\begin{center}
  \begin{minipage}[c]{0.5\textwidth}
    \includegraphics[width=\textwidth]{./results_reg/Square0_t1_8.eps.bz2}
  \end{minipage}%
  \begin{minipage}[c]{0.5\textwidth}
    \includegraphics[width=\textwidth]{./results_reg/Lshape0_t1_8.eps.bz2}
  \end{minipage}%
\end{center}

\begin{center}
  \begin{minipage}[c]{0.5\textwidth}
    \includegraphics[width=\textwidth]{./results_reg/Square1_t1_8.eps.bz2}
  \end{minipage}%
  \begin{minipage}[c]{0.5\textwidth}
    \includegraphics[width=\textwidth]{./results_reg/Lshape1_t1_8.eps.bz2}
  \end{minipage}%
\end{center}

\begin{center}
  \begin{minipage}[c]{0.5\textwidth}
    \includegraphics[width=\textwidth]{./results_reg/Square2_t1_8.eps.bz2}
  \end{minipage}%
  \begin{minipage}[c]{0.5\textwidth}
    \includegraphics[width=\textwidth]{./results_reg/Lshape2_t1_8.eps.bz2}
  \end{minipage}%
\end{center}
  \caption{Discrete solutions on different meshes at time $t=1.8$}
  \label{fig:ds}
\end{figure}

Except for different numerical dispersions, both edge elements and
vertex based elements largely produce the same solution on the
square. Conversely, there is a conspicuous qualitative difference of
the solutions in the presence of reentrant corners. A similar message
is sent by the temporal variation of electric and magnetic energies
plotted in figures \ref{fig:enreg0}, \ref{fig:enreg1}, \ref{fig:enreg2}. 

\begin{figure}[!htbp]
  \centering
  \begin{minipage}[c]{0.5\textwidth}
    \includegraphics[width=\textwidth]{./results_reg/Square0_en.eps.bz2}
  \end{minipage}%
  \begin{minipage}[c]{0.5\textwidth}
    \includegraphics[width=\textwidth]{./results_reg/Lshape0_en.eps.bz2}
  \end{minipage}%
  \caption{Energies, regularized: coarse
  meshes (left: square, right: L-shaped domain)}
  \label{fig:enreg0}
\end{figure}

\begin{figure}[!htbp]
  \centering
  \begin{minipage}[c]{0.5\textwidth}
    \includegraphics[width=\textwidth]{./results_reg/Square1_en.eps.bz2}
  \end{minipage}%
  \begin{minipage}[c]{0.5\textwidth}
    \includegraphics[width=\textwidth]{./results_reg/Lshape1_en.eps.bz2}
  \end{minipage}%
  \caption{Energies, regularized: intermediate meshes (left: square, right: L-shaped)}
  \label{fig:enreg1}
\end{figure}

\begin{figure}[!htbp]
  \centering
  \begin{minipage}[c]{0.5\textwidth}
    \includegraphics[width=\textwidth]{./results_reg/Square2_en.eps.bz2}
  \end{minipage}%
  \begin{minipage}[c]{0.5\textwidth}
    \includegraphics[width=\textwidth]{./results_reg/Lshape2_en.eps.bz2}
  \end{minipage}%
  \caption{Energies, regularized: fine
  meshes (left: square, right: L-shaped domain)}
  \label{fig:enreg2}
\end{figure}

\subsection{Discrete evolution without regularization}

Compared to the previous experiments all regularizing terms were
dropped. Snapshots of the discrete solutions are displayed for
$t=1.8$ in figures \ref{fig:dscurl}.

\begin{figure}[htbp]
\begin{center}
  \begin{minipage}[c]{0.5\textwidth}
    \includegraphics[width=\textwidth]{./results_curl/Square0_t1_8.eps.bz2}
  \end{minipage}%
  \begin{minipage}[c]{0.5\textwidth}
    \includegraphics[width=\textwidth]{./results_curl/Lshape0_t1_8.eps.bz2}
  \end{minipage}%
\end{center}

\begin{center}
  \begin{minipage}[c]{0.5\textwidth}
    \includegraphics[width=\textwidth]{./results_curl/Square1_t1_8.eps.bz2}
  \end{minipage}%
  \begin{minipage}[c]{0.5\textwidth}
    \includegraphics[width=\textwidth]{./results_curl/Lshape1_t1_8.eps.bz2}
  \end{minipage}%
\end{center}

\begin{center}
  \begin{minipage}[c]{0.5\textwidth}
    \includegraphics[width=\textwidth]{./results_curl/Square2_t1_8.eps.bz2}
  \end{minipage}%
  \begin{minipage}[c]{0.5\textwidth}
    \includegraphics[width=\textwidth]{./results_curl/Lshape2_t1_8.eps.bz2}
  \end{minipage}%
\end{center}
  \caption{Discrete solutions on different meshes at time $t=1.8$}
  \label{fig:dscurl}
\end{figure}

The overall appearence of the discrete solutions is strikingly similar
for both the edge element and nodal scheme regardless of the domain.
However, the nodal solution is marred by the presence of wiggles, which
are probably due to excited near-kernel modes. Nevertheless the energies
show good agreement up to $t=2.5$ until inevitable numerical dispersion
interferes.

\begin{figure}[!htbp]
  \centering
  \begin{minipage}[c]{0.5\textwidth}
    \includegraphics[width=\textwidth]{./results_curl/Square0_en.eps.bz2}
  \end{minipage}%
  \begin{minipage}[c]{0.5\textwidth}
    \includegraphics[width=\textwidth]{./results_curl/Lshape0_en.eps.bz2}
  \end{minipage}%
  \caption{Energies, plain scheme: coarse
  meshes (left: square, right: L-shaped domain)}
  \label{fig:encurl0}
\end{figure}

\begin{figure}[!htbp]
  \centering
  \begin{minipage}[c]{0.5\textwidth}
    \includegraphics[width=\textwidth]{./results_curl/Square1_en.eps.bz2}
  \end{minipage}%
  \begin{minipage}[c]{0.5\textwidth}
    \includegraphics[width=\textwidth]{./results_curl/Lshape1_en.eps.bz2}
  \end{minipage}%
  \caption{Energies, plain scheme: intermediate meshes (left: square, right: L-shaped )}
  \label{fig:encurl1}
\end{figure}

\begin{figure}[!htbp]
  \centering
  \begin{minipage}[c]{0.5\textwidth}
    \includegraphics[width=\textwidth]{./results_curl/Square2_en.eps.bz2}
  \end{minipage}%
  \begin{minipage}[c]{0.5\textwidth}
    \includegraphics[width=\textwidth]{./results_curl/Lshape2_en.eps.bz2}
  \end{minipage}%
  \caption{Energies, plain scheme: fine
  meshes (left: square, right: L-shaped domain)}
  \label{fig:encurl2}
\end{figure}

\bibliographystyle{siam}
\begin{footnotesize}
\bibliography{lit}
\end{footnotesize}

\end{document}