numerics.lyx

#LyX 2.3 created this file. For more info see http://www.lyx.org/
\lyxformat 544
\begin_document
\begin_header
\save_transient_properties true
\origin unavailable
\textclass article
\use_default_options true
\maintain_unincluded_children false
\language english
\language_package default
\inputencoding auto
\fontencoding global
\font_roman "default" "default"
\font_sans "default" "default"
\font_typewriter "default" "default"
\font_math "auto" "auto"
\font_default_family default
\use_non_tex_fonts false
\font_sc false
\font_osf false
\font_sf_scale 100 100
\font_tt_scale 100 100
\use_microtype false
\use_dash_ligatures true
\graphics default
\default_output_format default
\output_sync 0
\bibtex_command default
\index_command default
\paperfontsize default
\spacing single
\use_hyperref false
\papersize default
\use_geometry true
\use_package amsmath 1
\use_package amssymb 1
\use_package cancel 1
\use_package esint 1
\use_package mathdots 1
\use_package mathtools 1
\use_package mhchem 1
\use_package stackrel 1
\use_package stmaryrd 1
\use_package undertilde 1
\cite_engine basic
\cite_engine_type default
\biblio_style plain
\use_bibtopic false
\use_indices false
\paperorientation portrait
\suppress_date false
\justification true
\use_refstyle 1
\use_minted 0
\index Index
\shortcut idx
\color #008000
\end_index
\leftmargin 3cm
\topmargin 3cm
\rightmargin 3cm
\bottommargin 3cm
\secnumdepth 3
\tocdepth 3
\paragraph_separation indent
\paragraph_indentation default
\is_math_indent 0
\math_numbering_side default
\quotes_style english
\dynamic_quotes 0
\papercolumns 1
\papersides 1
\paperpagestyle default
\tracking_changes false
\output_changes false
\html_math_output 0
\html_css_as_file 0
\html_be_strict false
\end_header

\begin_body

\begin_layout Standard

\lang british
\begin_inset FormulaMacro
\newcommand{\tht}{\vartheta}
\end_inset


\begin_inset FormulaMacro
\newcommand{\ph}{\varphi}
\end_inset


\begin_inset FormulaMacro
\newcommand{\balpha}{\boldsymbol{\alpha}}
\end_inset


\begin_inset FormulaMacro
\newcommand{\btheta}{\boldsymbol{\theta}}
\end_inset


\begin_inset FormulaMacro
\newcommand{\bJ}{\boldsymbol{J}}
\end_inset


\begin_inset FormulaMacro
\newcommand{\bGamma}{\boldsymbol{\Gamma}}
\end_inset


\begin_inset FormulaMacro
\newcommand{\bOmega}{\boldsymbol{\Omega}}
\end_inset


\begin_inset FormulaMacro
\newcommand{\d}{\text{d}}
\end_inset


\begin_inset FormulaMacro
\newcommand{\t}[1]{\text{#1}}
\end_inset


\begin_inset FormulaMacro
\newcommand{\m}{\text{m}}
\end_inset


\begin_inset FormulaMacro
\newcommand{\v}[1]{\boldsymbol{#1}}
\end_inset


\end_layout

\begin_layout Standard

\lang british
\begin_inset FormulaMacro
\renewcommand{\t}[1]{\mathbf{#1}}
\end_inset


\end_layout

\begin_layout Section
Global interpolation techniques
\end_layout

\begin_layout Subsection
Series expansion in a Hilbert space
\end_layout

\begin_layout Standard
Consider a complete basis 
\begin_inset Formula $\{\ph_{n}\}$
\end_inset

 for functions in a Hilbert space with inner product 
\begin_inset Formula $\left\langle f,g\right\rangle $
\end_inset

.
 This means we can write every function as
\end_layout

\begin_layout Standard
\begin_inset Formula 
\begin{equation}
f(x)=\sum_{n=0}^{\infty}f_{n}\ph_{n}(x).
\end{equation}

\end_inset

To find the coefficients 
\begin_inset Formula $f_{n}$
\end_inset

 we do a projection
\begin_inset Formula 
\begin{equation}
\left\langle \ph_{m},f\right\rangle =\left\langle \ph_{m},\sum_{n=0}^{\infty}f_{n}\ph_{n}\right\rangle =\sum_{n=0}^{\infty}\left\langle \ph_{m},\ph_{n}\right\rangle f_{n}.
\end{equation}

\end_inset

Here we could move the infinite sum and the coefficients 
\begin_inset Formula $f_{n}$
\end_inset

 out of the inner product, as the latter is linear and continuous in both
 arguments.
 To solve this problem for 
\begin_inset Formula $f_{n}$
\end_inset

 we can truncate the sum with 
\begin_inset Formula $m,n\leq N$
\end_inset

 and solve a set of 
\begin_inset Formula $N+1$
\end_inset

 linear equations
\begin_inset Formula 
\begin{align}
\sum_{n=0}^{N}A_{mn}f_{n} & =b_{m},\quad\text{where}\\
A_{mn} & =\left\langle \ph_{m},\ph_{n}\right\rangle ,\quad b_{m}=\left\langle \ph_{m},f\right\rangle .
\end{align}

\end_inset

If the property of 
\series bold
\emph on
orthogonality 
\end_layout

\begin_layout Standard
\begin_inset Formula 
\[
\left\langle \ph_{m},\ph_{n}\right\rangle =\left\langle \ph_{n},\ph_{n}\right\rangle \delta_{mn}
\]

\end_inset

is fulfilled we don't need to limit ourself to a finite system, but can
 solve exactly
\begin_inset Formula 
\begin{equation}
f_{n}=\frac{\left\langle \ph_{n},f\right\rangle }{\left\langle \ph_{n},\ph_{n}\right\rangle }.\label{eq:coeff}
\end{equation}

\end_inset

Two bases that are orthogonal under the usual 
\begin_inset Formula $L_{2}$
\end_inset

 norm
\begin_inset Formula 
\begin{equation}
\left\langle f,g\right\rangle =\int_{a}^{b}f(x)g(x)\,\d x
\end{equation}

\end_inset

over some finite interval 
\begin_inset Formula $[a,b]$
\end_inset

 are Fourier (sine/cosine/complex exponential) series bases and Legendre
 polynomials.
 There exists a wider class of 
\emph on
orthogonal polynomials 
\emph default
that are orthogonal with respect to an inner product with a weight function
 
\begin_inset Formula $w$
\end_inset

 with
\begin_inset Formula 
\begin{equation}
\left\langle f,g\right\rangle =\int_{a}^{b}w(x)f(x)g(x)\,\d x,
\end{equation}

\end_inset

where either 
\begin_inset Formula $a$
\end_inset

 or 
\begin_inset Formula $b$
\end_inset

 or both limits can be infinite.
\end_layout

\begin_layout Subsection
Fourier series
\end_layout

\begin_layout Standard
TODO
\end_layout

\begin_layout Subsection
Legendre polynomials
\end_layout

\begin_layout Standard
Every square-integrable function 
\begin_inset Formula $f$
\end_inset

 on the interval 
\begin_inset Formula $[0,1]$
\end_inset

 can be represented as
\begin_inset Formula 
\begin{equation}
f(x)=\sum_{n=0}^{\infty}f_{n}P_{n}(x).
\end{equation}

\end_inset

The coefficients 
\begin_inset Formula $f_{n}$
\end_inset

 can be computed by using the especially simple orthogonality relation with
\begin_inset Formula 
\begin{equation}
\left\langle P_{m},P_{n}\right\rangle =\int_{-1}^{1}P_{m}(x)P_{n}(x)\d x=\frac{2}{2n+1}\delta_{mn}.
\end{equation}

\end_inset

According to Eq.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand eqref
reference "eq:coeff"

\end_inset

 coefficients in the Legendre basis are given by
\begin_inset Formula 
\begin{equation}
f_{n}=\frac{2n+1}{2}\int_{-1}^{1}P_{n}(x)f(x)\d x.
\end{equation}

\end_inset


\end_layout

\begin_layout Subsubsection
Derivatives
\end_layout

\begin_layout Standard
We can compute derivatives of Legendre polynomials by
\begin_inset Formula 
\begin{equation}
(n-m)P_{n}^{(m)}(x)=(2n-1)xP_{n}^{(m)}(x)-(n-1+m)P_{n-2}^{(m)}(x)
\end{equation}

\end_inset

and
\begin_inset Formula 
\begin{align}
P_{m}^{(m)}(x) & =\frac{(2m)!}{2^{m}m!},\\
P_{k}^{(m)}(x) & =0\quad\text{for }k<m.
\end{align}

\end_inset

In particular, first derivatives are given by
\begin_inset Formula 
\begin{equation}
(n-1)P_{n}^{\prime}(x)=(2n-1)xP_{n}^{\prime}(x)-(n-1+m)P_{n-2}^{\prime}(x)
\end{equation}

\end_inset

starting the recursion with
\begin_inset Formula 
\begin{align}
P_{0}^{\prime}(x) & =0,\\
P_{1}^{\prime}(x) & =1.
\end{align}

\end_inset


\end_layout

\begin_layout Subsection
Hierarchical combination of different bases
\end_layout

\begin_layout Standard
The described method seems to be closely related to hierarchical and multigrid/m
ultilevel methods, and coarse grid preconditioners.
 The advantage of this method is that we can 
\begin_inset Quotes eld
\end_inset

filter out
\begin_inset Quotes erd
\end_inset

 coarse scales of the problem by a global basis and treat the remaining
 fine scales in a local basis.
 This way we get the 
\begin_inset Quotes eld
\end_inset

best of both worlds
\begin_inset Quotes erd
\end_inset

 and retain global continuity.
\end_layout

\begin_layout Standard
Say we have two (incomplete/finite) bases 
\begin_inset Formula $\{\ph_{n}^{(1)}\}$
\end_inset

 and 
\begin_inset Formula $\{\ph_{k}^{(2)}\}$
\end_inset

.
 We want to solve a linear differential equations
\begin_inset Formula 
\begin{equation}
Df(x)=g(x)\label{eq:D}
\end{equation}

\end_inset

on some domain.
 Now we first solve a discrete equation in the coefficients with
\begin_inset Formula 
\begin{equation}
\sum_{n}D_{mn}^{(1)}f_{n}^{(1)}=g_{m}^{(1)},\label{eq:D1}
\end{equation}

\end_inset

where
\begin_inset Formula 
\begin{equation}
f^{(1)}(x)=\sum_{n}f_{n}^{(1)}\ph_{n}(x),\quad g^{(1)}(x)=\sum_{n}g_{n}^{(1)}\ph_{n}(x).
\end{equation}

\end_inset

Solving Eq.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand eqref
reference "eq:D1"

\end_inset

 leads an approximate solution of Eq.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand eqref
reference "eq:D"

\end_inset

 but yields the exact relation
\begin_inset Formula 
\begin{equation}
Df^{(1)}(x)=g^{(1)}(x).\label{eq:D-1}
\end{equation}

\end_inset

In the second step we subtract the solution from the right-hand side and
 try to solve
\begin_inset Formula 
\begin{equation}
Df(x)=g(x)-g^{(1)}(x)\equiv g^{(2)}(x).
\end{equation}

\end_inset

Thus
\begin_inset Formula 
\begin{equation}
\sum_{n}D_{jk}^{(2)}f_{k}^{(2)}=g_{j}^{(2)}.
\end{equation}

\end_inset

Again the solution is not exact but we obtain rather
\begin_inset Formula 
\[
Df^{(2)}(x)=g^{(2)}(x)-g^{(\text{res})}(x)
\]

\end_inset

with some residual 
\begin_inset Formula $g^{(\text{res})}$
\end_inset

.
\end_layout

\begin_layout Standard
To sum up, we have split the problem into 
\begin_inset Formula 
\begin{align}
Df^{(1)}(x) & =g^{(1)}(x)\equiv g(x)-g^{(2)}(x),\\
Df^{(2)}(x) & =g^{(2)}(x)-g^{(\text{res})}(x),
\end{align}

\end_inset

where the first equation is exactly solved since 
\begin_inset Formula $g^{(2)}$
\end_inset

 is set to the residual with respect to the exact solution.
 The approximation appears at the very bottom with 
\begin_inset Formula $g^{(\text{res})}(x)$
\end_inset

.
 If we take a sum we see that 
\begin_inset Formula $g^{(2)}$
\end_inset

 cancels and the overall approximate solution
\begin_inset Formula 
\[
f^{(h)}=f^{(1)}+f^{(2)}
\]

\end_inset

fulfills
\begin_inset Formula 
\begin{align}
Df^{(h)}(x) & =g(x)-g^{(\text{res})}(x).
\end{align}

\end_inset

We could continue this process for an arbitrary number of hierarchy levels
 by setting 
\begin_inset Formula $g^{(\text{res})}$
\end_inset

 to 
\begin_inset Formula $g^{(3)}$
\end_inset

 and so on.
\end_layout

\begin_layout Section
Root finding
\end_layout

\begin_layout Subsection
Newton's method
\end_layout

\begin_layout Standard
When sufficiently close to a local minimum we use a Taylor expansion
\begin_inset Formula 
\begin{equation}
0=f(x_{0})=f(x)+(x_{0}-x)f^{\prime}(x)+\mathcal{O}(|x-x_{0}|^{2})
\end{equation}

\end_inset

to obtain the next best guess 
\begin_inset Formula $x_{k+1}\approx x_{0}$
\end_inset

 via
\begin_inset Formula 
\begin{equation}
x_{k+1}=x_{k}-\frac{f(x_{k})}{f^{\prime}(x_{k})}.\label{eq:xkp1}
\end{equation}

\end_inset

Or in 
\begin_inset Formula $n$
\end_inset

 dimensions
\begin_inset Formula 
\begin{equation}
0=\boldsymbol{F}(\boldsymbol{x}_{0})=\boldsymbol{F}(\boldsymbol{x})+(D_{\boldsymbol{x}}\boldsymbol{F}(\boldsymbol{x}))(\boldsymbol{x}-\boldsymbol{x}_{0})+\mathcal{O}(|\boldsymbol{x}-\boldsymbol{x}_{0}|^{2})
\end{equation}

\end_inset

The equivalent iteration rule to Eq.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand eqref
reference "eq:xkp1"

\end_inset

 would be
\begin_inset Formula 
\begin{equation}
\boldsymbol{x}_{k+1}=\boldsymbol{x}_{k}-(D_{\boldsymbol{x}}\boldsymbol{F}(\boldsymbol{x}))^{-1}\boldsymbol{F}(\boldsymbol{x}_{k}).
\end{equation}

\end_inset

Instead of inverting the Jacobian 
\begin_inset Formula $D_{\boldsymbol{x}}\boldsymbol{F}$
\end_inset

 here it is often more efficient and numerically stable to solve the linear
 system
\begin_inset Formula 
\begin{equation}
(D_{\boldsymbol{x}}\boldsymbol{F}(\boldsymbol{x}))\Delta\boldsymbol{x}_{k+1}=-\boldsymbol{F}(\boldsymbol{x}_{k})
\end{equation}

\end_inset

in 
\begin_inset Formula $\Delta\boldsymbol{x}_{k+1}\equiv\boldsymbol{x}_{k+1}-\boldsymbol{x}_{k}$
\end_inset

 instead.
\end_layout

\begin_layout Standard
Modified versions of Newton's method are
\end_layout

\begin_layout Itemize
Inexact Newton: for inverting the Jacobian, use an interative method (e.g.
 GMRES) and stop before convergence reached.
\end_layout

\begin_layout Itemize
Quasi-Newton: use an approximate Jacobian instead of the exact one (Broyden,
 Powell hybrid, etc.)
\end_layout

\begin_layout Subsection
Linear extrapolation
\end_layout

\begin_layout Standard
Since 
\begin_inset Formula $\boldsymbol{F}$
\end_inset

 and 
\begin_inset Formula $\nabla\boldsymbol{F}$
\end_inset

 are only evaluated at the second-last position 
\begin_inset Formula $\boldsymbol{x}_{N-1}$
\end_inset

 before convergence is reached at 
\begin_inset Formula $\boldsymbol{x}_{N}$
\end_inset

, we cannot directly re-use quantities that have been computed during this
 evaluation without another evaluation at 
\begin_inset Formula $\boldsymbol{x}_{N}$
\end_inset

.
 Consider a quantity 
\begin_inset Formula $a(\boldsymbol{x})$
\end_inset

.
 In case both 
\begin_inset Formula $a(\boldsymbol{x}_{N-1})$
\end_inset

 and 
\begin_inset Formula $\nabla a(\boldsymbol{x}_{N-1})$
\end_inset

 are known from the last iteration step, a good approximation at 
\begin_inset Formula $\boldsymbol{x}_{N}$
\end_inset

 is given by the linear extrapolation
\begin_inset Formula 
\begin{equation}
a(\boldsymbol{x}_{N})=a(\boldsymbol{x}_{N-1})+(\boldsymbol{x}_{N}-\boldsymbol{x}_{N-1})\cdot\nabla a(\boldsymbol{x}_{N-1})+\mathcal{O}(|\boldsymbol{x}_{N}-\boldsymbol{x}_{N-1}|^{2}).
\end{equation}

\end_inset

Thus we can save one evaluation of 
\begin_inset Formula $a$
\end_inset

 in computations done at the final position 
\begin_inset Formula $\boldsymbol{x}_{N}$
\end_inset

 after convergence of a root-finding scheme.
\end_layout

\begin_layout Subsection
Error analysis
\end_layout

\begin_layout Standard
See also http://www.math.usm.edu/lambers/mat460/fall09/lecture12.pdf and http://www.m
ath.niu.edu/~dattab/MATH435.2013/ROOT_FINDING.pdf
\end_layout

\begin_layout Standard
Suppose we are performing a fixed-point iteration
\begin_inset Formula 
\begin{equation}
x_{k+1}=g(x_{k})
\end{equation}

\end_inset

converging towards 
\begin_inset Formula $x^{\star}$
\end_inset

 with 
\begin_inset Formula $f(x^{\star})=0$
\end_inset

.
 In practice we stop either if either
\begin_inset Formula 
\begin{equation}
|f(x_{n})|<\text{atol}
\end{equation}

\end_inset

or
\begin_inset Formula 
\begin{equation}
\frac{|x_{n}-x_{n-1}|}{|x_{n}|}<\text{rtol}
\end{equation}

\end_inset

or the number of iteration has reached its maximum.
\end_layout

\begin_layout Section
Numerical solution of differential equations
\end_layout

\begin_layout Itemize
Particle methods move particles along characteristics of PDE and/or update
 weights
\end_layout

\begin_layout Itemize
FEM 
\begin_inset Quotes eld
\end_inset

particles
\begin_inset Quotes erd
\end_inset

 live on local support at fixed positions
\end_layout

\begin_layout Itemize
Kernel-based interpolation / smoothing / meshless methods are a 
\begin_inset Quotes eld
\end_inset

hybrid
\begin_inset Quotes erd
\end_inset

 of particle and FEM
\end_layout

\begin_layout Itemize
FVM
\end_layout

\begin_layout Itemize
DG
\end_layout

\begin_layout Itemize
time-stepping for ODEs
\end_layout

\end_body
\end_document