Numerical Methods
for PDE

General facts about PDEs

Examples of important elementary PDEs

Existence and uniqueness of solutions

This theorem though fails if any of the function is just simply smooth rather than analytic!
see https://en.wikipedia.org/wiki/Lewy's_example.

Classification

Basic examples:

Classification

• Elliptic PDEs describe processes that have already reached a steady state, and so are time-independent.
• Hyperbolic PDEs describe time-dependent conservative processes, such as convection, thsat are not evolving towards a steady state.
• Parabolic PDEs describe time-dependent dissipative processes, such as diffusion, that are evolving towards a steady state.

Hyperbolic PDEs

Advection equation

Continuity equation

Advection equation

1D Advection equation

Characteristics

Solving Advection PDE with FDM

$$t_n = n\Delta t,\;0\leq n\leq N;\;\;x_m = m\Delta x,\;0\leq m\leq M$$

and $u_t$ and $u_x$ with one of their discrete approximations.

1. FTBS (Upwind Method)

FTBS Diagram

FTBS Initial Conditions

FTBS Initial Conditions

Implementation in Python

Method convergence: The CFL condition

The CFL condition

The CFL condition

Method convergence: The von Neumann condition

2. FTCS

FTCS Diagram

FTCS vc FTBS

3. The Lax Method

Lax Diagram

Lax's CFL and Von Neumann conditions

General properties of numerical pde methods

Consistency, stability, convergence

• Consistent if its truncation error (i.e. the local error) goes to zero as $\Delta t,\Delta x\to0$.
• Stable if the operator $T$ sending initial conditions to the corresponding numerical solutions is bounded, namely there are no initial conditions on which the numerical solution diverges (von Neumann analysis).
• Convergent if the global error goes to zero as $\Delta t,\Delta x\to0$.

The Lax Equivalence Theorem

Parabolic PDEs

Diffusion equation

Heat equation

1D Heat equation

1. FTCS

FTCS Stability

FTCS accuracy problem

Implementation in MATLAB

2. CTCS (Richardson Method)

3. DuFort-Frankel method

3. DuFort-Frankel method

3. DuFort-Frankel method

4. BTCS method

$$ u_{n,j} = u_{n-1,j} + \nu\left(u_{n,j+1}-2u_{n,j}+u_{n,j-1}\right)\,, $$

that can be solved "backwards":

$$ u_{n-1,j} = u_{n,j} - \nu\left(u_{n,j+1}-2u_{n,j}+u_{n,j-1}\right)\,, $$

4. BTCS Solution

4. BTCS - Stability

so that $$ G(k) = \frac{1}{1+4\nu\sin^2\frac{k\Delta x}{2}} $$

This method therefore is stable for all $\Delta t$ and $\Delta x$!

4. BTCS - Weak point

Implementation in MATLAB

5. Crank-Nicholson method

5. Crank-Nicholson Stability

$$ G(k) = \frac{1-2r\sin^2\frac{k\Delta x}{2}}{1+2r\sin^2\frac{k\Delta x}{2}} $$

This method therefore is also unconditionally stable
plus
is 2nd order in both time and space.

Implementation in MATLAB

Elliptic PDEs

Poisson Equation

$$ u_{xx} + u_{yy} = \rho\;\;\hbox{ in }S\,,\;\;u|_{\partial S}=f, $$

is named Poisson Equation or, when $f=0$, Laplace Equation.

Example: electrostatic electric potential

Electric Potential of a point-like charge

Finite Differences -- Explicit method

Finite Differences -- Explicit method

We set $\Delta x=\Delta y=h$, so we get the following system:

$$ V_{i+1,j}+V_{i,j+1}-4V_{i,j}+V_{i-1,j}+V_{i,j-1}=4\pi h^2\rho_{i,j} $$

This Finite Differences scheme can be summarized
by the 5-points stencil diagram below:

Finite Differences -- Explicit method

$$ M = \begin{pmatrix} B&-I&&&&&0\\ -I& B&-I&&&&\\ &.&.&.&&&\\ &&.&.&.&&\\ &&&.&.&.&\\ &&&&-I&B&-I\\ 0&&&&&-I&B\\ \end{pmatrix} $$

is a $N^2$ block matrix, where...

Finite Differences -- Explicit method

Local Truncation Error

Global Error

Stability and Accuracy

Stability and Accuracy

Stability and Accuracy

Evaluating in concrete a numerical solution

Solving $M\vec V=\vec R$ -- LU method

The complexity of this method is $O\left(N^3\right)$. For banded matrices of bandwidth $d$, the complexity goes down to $O\left(d^2N\right)$.

Solving $M\vec V=\vec R$ -- Jacobi method

First Implementation in Python

Second Implementation in Python

Solving $M\vec V=\vec R$ -- Gauss-Seidel method

Solving $M\vec V=\vec R$ -- GS with SOR relaxation

The idea behind the

Finite Elements Method

A Linear system

We want to solve the equation $$ A v = r\,. $$

Recall that every finite dimensional vector space admits bases $\{e_1,\dots,e_n\}$, $n=\dim V$,
so that every $v\in V$ can be written in a unique way as $$v=\sum_{i=1}^n v_i e_i\,.$$

An equivalent formulation

$\langle A v,e_i\rangle = \langle r,e_i\rangle\,$, for all $1\leq i\leq n$.

An equivalent formulation

In turn, this is the same as asking

$\langle A v,s\rangle = \langle r,s\rangle\,$, for all $s\in V$.

While in finite dimension these two formulations are equivalent,
in infinite dimension the latter is weaker,
namely it has in principle more solutions than the original one.

Example 1 -- ODEs

The completion of $V$ with respect to $\langle,\rangle$ is called $L^2[0,1]$ and is the space of all measurable functions $u$ such that $$\|u\|^2=\int_0^1|u(x)|^2 dx<\infty\,.$$

$$ -u''(x)+r(x)u(x)=f(x)\,,\;x\in[0,1]\,,\;\;u(0)=u(1)=0\,, $$

with $f,r\in V$. This is a linear equation in $V$ and can be written as $$ Au=f\,,\;u\in C^\infty_0[0,1] $$ where $A$ is the linear operator from $V$ to itself $$ -\frac{d^2}{dx^2}+r(x) $$ and $C^\infty_0[0,1]$ is the subset of smooth functions that are zero at $x=0,1$.

In this equation all functions appear under integral sign,
so there is no actual need to require so much regularity:

it is enough to ask that both $u$ and $s$ belong to $H^1_0[0,1]$,
the space of square-integrable functions whose first derivative
is also square-integrable and that converge to $0$ for $x\to0^+,1^-$.

Approximating weak solutions

If $\{\phi_i\}$ is a basis for the subspace $H^1_0[0,1]$ of all functions of $L^2[0,1]$ with square-integrble gradient and which tend zero for $x\to0^+$ and $x\to1^-$, then

$\langle u',\phi'_i\rangle+r\langle u,\phi_i\rangle=\langle f,\phi_i\rangle$

for all $i$. Since $u=\sum_{i=1}^\infty u_i\phi_i$ and $f=\sum_{i=1}^\infty f_i\phi_i$ we have that

$\sum_{j=1}^\infty u_j\left( \langle \phi'_j,\phi'_i\rangle + r\langle \phi_j,\phi_i\rangle\right) = \sum_{j=1}^\infty f_j\langle \phi_j,\phi_i\rangle$

By truncating this sum to some $N$, we get that $$ AU=F $$ form some matrix $A$ and $F$.

Example 2 -- PDEs

Again, $V$ is an infinite dimensional vector space with scalar product

$$ \langle f,g \rangle = \displaystyle\iint_S\!\!\!f(x,y)g(x,y)dxdy\,. $$

$$ -u_{xx}(x,y)-u_{yy}(x,y)=f(x,y)\,,\;x\in[0,1]\,,\;\;u=0 on \partial S\,, $$ where $f\in V$ and $\partial S$ is the boundary of $S$.

This is a linear equation in $V$ and can be written as $$ Au=f\,,\;u\in C^\infty_0(S) $$ where $A$ is the linear operator from $V$ to itself $$ -(\partial_{xx}+\partial_{yy}) $$ and $C^\infty_0(S)$ is the subset of smooth functions that are zero at the boundary of $S$.

Then we integrate by parts the lhs terms as above: $$ \langle u_{xx},s\rangle = -\langle u_x,s_x\rangle\,,\;\langle u_{yy},s\rangle = -\langle u_y,s_y\rangle $$ The equation now is

$$ \langle\nabla u,\nabla s\rangle = \langle f,s\rangle\,,\;\hbox{ for all }s\in C^\infty_0(S) $$

Again, in this equation all functions appear under integral sign so there is no actual need to require so much regularity: it is enough to ask that both $u$ and $s$ belong to $H^1_0(S)$, the space of $L^2$ functions on $D$ whose gradient is also $L^2$ and that converge to 0 on any sequence of points converging to the boundary.

Approximating weak solutions

If $\{\phi_i\}$ is a basis for $H^1_0(S)$, then

$\langle \nabla u,\nabla \phi_i\rangle=\langle f,\phi_i\rangle$

for all $i$. Since $u=\sum_{i=1}^\infty u_i\phi_i$ and $f=\sum_{i=1}^\infty f_i\phi_i$ we have that

$\sum_{j=1}^\infty u_j \langle \nabla \phi_j,\nabla\phi_i\rangle = \sum_{j=1}^\infty f_j\langle\phi_j,\phi_i\rangle$

By truncating this sum to some $N$, we get that $$ AU=F $$ form some matrix $A$ and $F$.

Implementation of the

Finite Elements Method

Why FEM?

This is accomplished mainly with the following two points:

1. Generating a suitable mesh (possibly tailored on the specifics of the PDE) of $D$ with $N$ nodes $\{p_1,\dots,p_N\}$ and a corresponding set of functions $\{\phi_1,\dots,\phi_N\}$ such that:
   1. $\phi_i(p_j)=\delta_{i,j}$;
   2. $\phi_i(p_i)=0$ on every mesh polygon not containing $p_i$.
   3. $\phi_i$ has an absolute max at $p_i$ and no other local max or saddle point on the polygons on which is not zero.
2. Using the PDE's weak version and restricting the eq. to $\hbox{span}\,\{\phi_1,\dots,\phi_N\}\subset H^1(D)$.

Meshes

1. The union $M_D$ of all $P_i$'s is an approximation of $D$, namely the measure of the set of points not in $M_D\cap D$ is small with respect to the measure of $D$.
2. For any $i\neq j$, $P_i$ and $P_j$ are either disjoint or share a whole side; in particular, no vertex of one polygon falls in the interior of the side of another one.

Mesh generation

			      
 m = 10;
 n = 10;
 [x,y] = ndgrid((0:m-1)/(m-1),(0:n-1)/(n-1));
 pts = [x(:),y(:)];
 tri = [1,2,m+2; 1,m+2,m+1];
 tri = kron(tri,ones(m-1,1))+kron(ones(size(tri)),(0:m-2)');
 tri = kron(tri,ones(n-1,1))+kron(ones(size(tri)),(0:n-2)'*m);
 trimesh(t,pts(:,1),pts(:,2),zeros(N,1))
			      
			      

Functions Bases

In this type of mesh, a reasonable set of $\phi_i$ is the set of piece-wise linear functions that, on any triangle $T$ of the mesh, are either identically zero or linear.

Let $(x_i,y_i)$, $i=1,2,3$ be the vertices of $T$ and $\psi_i=a_i+b_ix+c_iy$, $i=1,2,3$, the restriction to $T$ of the three basis functions that are not identically zero on $T$ and such that $\psi_i(x_j,y_j)=\delta_{ij}$.

Functions Bases

namely

$$ \begin{pmatrix} a_1&a_2&a_3\cr b_1&b_2&b_3\cr c_1&c_2&c_3\cr \end{pmatrix} = M^{-1}_T $$

Implementation in MATLAB

Spectral methods

Fourier transform

Using Fourier Transform to solve ODEs

For example, consider $-u''(x)+u(x)=f(x)$. In momentum space the ODE just becomes the algebraic equation

$$ k^2\hat u(k)+\hat u(k)=\hat f(k) $$

so that $\hat u(k)=\hat f(k)/(1+k^2)$ and finally we get the ODE solution $$ u(x) = \frac{1}{2\pi}\int\limits_{-\infty}^\infty\frac{\hat f(k)}{1+k^2}\,e^{ikx}\,dk\,. $$

Fourier series

$$ u(x) = \sum_{k=-\infty}^\infty \hat u_k e^{2\pi ikx/L} $$

and

$$ \hat u_k = \frac{1}{L}\int\limits_{0}^L u(x) e^{-2\pi ikx/L}dx $$ so that $(u')_k = ik u_k$.

Clearly any DE under Fourier transform will become an algebraic equation, we can solve it and then use the inverse transformation to get the solution in $x$ after cutting the infinite series to $N$.

Discrete Fourier trasform

Suppose we want to solve a Dirichlet problem on $[0,L]$, so that $$ u(x)=\sum_{k=1}^\infty a_k \sin(\pi kx/L) $$

Then (Discrete Fourier Tranform) $$ u_l=u(l h)=\sum_{k=1}^N a_k \sin(\pi kl/(N+1))\,,\;1\leq l\leq N $$

is a $N\times N$ system in $a_1,\dots,a_N$.

Its solution is the Inverse Discrete Fourier Tranform $$ a_k=\frac{2}{N+1}\sum_{l=1}^N u_l\sin(\pi kl/(N+1))\,. $$

Von Neumann

Stability

Hence $|G(z)|>1$ for any value of $r$, namely this numerical method is unstable for any choice of $\Delta t$ and $\Delta x$!

A remedy: the Lax–Friedrichs method

$$ \frac{u_{n,j+1}+u_{n,j-1}}{2}\,, $$

so that the system becomes

$$ u_{n+1,j} = \frac{u_{n,j+1}+u_{n,j-1}}{2} + \frac{r}{2}\left(u_{n,j+1}-u_{n,j-1}\right)\,, $$

Under the guess $u(t,x) = e^{ik(x-ct)}$ we find that

$$ G(k) = \cos\left({k\Delta x}\right) + i\,r\cos\left({k\Delta x}\right) $$

so that $|G(k)|\leq1$ when $c\Delta t\leq\Delta x$.