DDG2016

Wave equation

Let $M$ be a surface and $f:\mathbb{R}\times M\to \mathbb{R}$ a function. The wave equation is given by:
$$\ddot{f}=\mathrm{\Delta }f.$$
After the space discretization we have $f:\mathbb{R}\times V\to \mathbb{R},$ $p:=\left(\begin{array}{c}{f}_1\\ ⋮\\ f_n\end{array}\right)\in {\mathbb{R}}^n$ and $\mathrm{\Delta }\in {\mathbb{R}}^{n\times n}$, where $n:=|V|$ denotes the number of vertices. Therefore, the wave equation becomes a linear system of second order ODE’s. In the last lecture we saw that the Verlet method works well for second order ODE’s. Using this for the Wave equation we get:$$\frac{p_{k-1}-2p_k+p_{k+1}}{{\delta }^2}=\mathrm{\Delta }{p}_k,$$
where $\delta$ is the time step.
For initial data $p_0,p_1$, the hope is that the energy of the system neither grows exponentially nor that the system comes to a halt. A classical example for the $1-$dimensional wave equation is an oscillating string. This leads in analogy to music instruments to two kinds of initial data:

1. Cembalo initial data: $$p_0=p_1,$$ here the string gets plucked, i.e. at time $t=0$ the oscillation is $p_o$ and the speed is zero.
2. Piano initial data: $$p_0=0,$$ here the string gets struck i.e.  at time $t=0$ the oscillation is zero but the speed is not.

A different way to find solutions of the wave equations is to calculate the eigenfunctions of the laplace operator, $g:M\to \mathbb{R}$:
$$\mathrm{\Delta }g=-{\omega }^{2}g.$$
Since the laplace operator is symmetric and negative definite, $\mathrm{\Delta }$ is diagonalizable and there exist eigenvalues ${\lambda }_i:=–{\omega }_i^2$ with corresponding eigenfunctions $g_i$ such that such that:
$$0=w_o^2<w_1^2<w_2^2\dots .$$
Note that $0$ zero is an eigenvalue, because the kernel of $\mathrm{\Delta }$ consists of the constant functions. If we computed the eigenfunctions, solutions of the wave equation can be constructed in the following way:
$$f(t,p):=\sum _{i=1}^n(a_i\cos ({\omega }_{i}t)+b_i\sin ({\omega }_{i}t))g_i(p),a_i,b_i\in \mathbb{R}.$$
The coefficients $a_i$, $b_i$ have to be chosen in a way such that the initial conditions are satisfied. Note that for different surfaces the spectrum of eigenvalues is different. A sphere “sounds” different then the bunny.  The disadvantage of this method is, that finding the eigenfunctions of the surface and computing the coefficients for the initial conditions can be very expensive, especially if the surface has many vertices.

Heat equation

For $f:\mathbb{R}\times M\to \mathbb{R}$, the heat equation is given by:
$$\dot{f}=\mathrm{\Delta }f.$$
We will start our investigation by considering the case $M={\mathbb{S}}^1.$ There we can describe a function $\tilde{f}:{\mathbb{S}}^1\to \mathbb{R}$ by an periodic function $f:\mathbb{R}\to \mathbb{R},$ with  $f(x+2\pi )=f(x).$

Our aim is it to use the Euler forward method to solve the heat equation on ${\mathbb{S}}^1$, therefore we discretize the space, i.e. the interval $[0,2\pi ]$, as follows:
$$x_k=\frac{2\pi k}{n},x_{k+1}=x_k+h,\text{with }h=\frac{2\pi }{n}.$$
In an earlier lecture we discussed discrete periodic functions and defined a canonical basis for them, such that we can write the components of $p_j=\left(\begin{array}{c}{u}_1\\ ⋮\\ u_{2n-1}\end{array}\right)$ as
$$u_k=\frac{1}{\sqrt{2\pi }}\sum _{m=-n}^n{a}_m{e}^{\frac{i2\pi km}{n}}{a}_m\in \mathbb{C}.$$
Since we are considering real valued functions there has to hold $a_m={\overline{a}}_{-m}$ for all $m\in \{-n,\dots ,0,\dots ,n\}$.
If we now use the Euler forward method we get:
$$p_{j+1}=p_j+\delta \mathrm{\Delta }{p}_k,$$
where $\delta$ denotes the time step size. The discrete laplacian on ${\mathbb{S}}^1$ can be considered as an approximation of the second spatial derivative and is defined as:
$${\left(\mathrm{\Delta }u\right)}_k=\frac{u_{k-1}-2u_k+u_{k+1}}{h^2}.$$
By linearity it is enough to consider the case $u_k=e^{\frac{i2\pi km}{n}}$:

$$\begin{array}{rl}{\tilde{u}}_k& =u_k+{\left(\mathrm{\Delta }u\right)}_k\\ & =e^{\frac{i2\pi km}{n}}+\delta \frac{e^{\frac{i2\pi (k-1)m}{n}}-2e^{\frac{i2\pi km}{n}}+e^{\frac{i2\pi (k+1)m}{n}}}{h^2}\\ & =(1+\delta \frac{e^{\frac{i2\pi m}{n}}+e^{\frac{-i2\pi m}{n}}+2}{h^2})e^{\frac{i2\pi km}{n}}\\ & =(1+\frac{2\delta }{h^2}(\cos \left(\frac{2\pi m}{n}\right)-1))u_k.\end{array}$$
The system will blow up after a certain time if there is an $m$ such that:
$$\begin{array}{rl}& |1+\frac{2\delta }{h^2}(\cos \left(\frac{2\pi m}{n}\right)-1)|>1\\ \iff & \frac{2\delta }{h^2}(\cos \left(\frac{2\pi m}{n}\right)-1)<-2.\end{array}$$
The left hand side reaches its minimum if $m=\frac{n}{2}$ i.e $\cos \left(\frac{2\pi m}{n}\right)=-1$ and we get:
$$\begin{array}{rl}& \frac{2\delta }{h^2}>1\\ \iff & \delta >\frac{h^2}{2}.\end{array}$$

This gives us a connection between the space and time discretization that is not wanted. If we use the Euler forward method to solve the heat equation on ${\mathbb{S}}^1$ and we refine the time step by the factor of $10$, we have to refine the space discretization by the factor of $50$, to avoid the solution to blow up. Therefore it makes more sense to use the Euler backward method for this problem.

Let $\mathrm{\Omega }\subset {\mathbb{R}}^n$ be a domain, we will consider $\mathrm{\Omega }$ to be the “space” and functions:
$$\begin{array}{r}f:\mathbb{R}\times \mathrm{\Omega }\to \mathbb{R},\\ (t,x)\mapsto f(t,x),\end{array}$$
will be view as time dependent functions
$$\begin{array}{rl}& f_t:\mathrm{\Omega }\to \mathbb{R},\\ & f_t(x):=f(t,x).\end{array}$$

We want to investigate partial differential equations  (PDE’s) for $f$, that involve the Laplace operator. Here the Laplace operator is defined just for the space coordinates $x$:
$$\left(\mathrm{\Delta }f\right)(t,x):=\left(\mathrm{\Delta }{f}_t\right)(x):=\sum _{i=1}^n\frac{{\partial }^2{f}_t(x)}{\partial x_i^2}.$$
Typical examples for such PDE’s from physics are:

1. Wave equation: $$\frac{{\partial }^{2}f}{\partial t^2}=c^2\mathrm{\Delta }f,c\in \mathbb{R}.$$
2. Heat equations: $$\frac{\partial f}{\partial t}=c\mathrm{\Delta }f,c\in \mathbb{R}.$$

We will investigate these PDS’s on closed triangulated surfaces  $\mathrm{\Omega }=M=(V,E,F)$, where we only have finitely many function values i.e. :
$$\begin{array}{rl}& f:\mathbb{R}\to {\mathbb{R}}^V,\\ & f(t)=\left(\begin{array}{c}{f}_1(t)\\ ⋮\\ f_N(t)\end{array}\right).\end{array}$$
On $M$ we have the usual Laplace operator $\mathrm{\Delta }:{\mathbb{R}}^V\to {\mathbb{R}}^V$ and for the above PDE’s we obtain ordinary differential equations (ODE’s):

1. Wave equation: $$\ddot{f}(t)=c^2\mathrm{\Delta }f(t).$$
2.  Heat equations: $$\dot{f}(t)=c\mathrm{\Delta }f(t).$$

This is a big improvement since it is fare more complicated to solve PDE’s then ODE’s.
Before we will present some methods to solve ODE’S numerically consider the Pendulum equation as example:
$$\ddot{\phi }=-\sin (\phi ).$$

The standard procedure to solve a second order ODE is to transfer it in a system of two ODE’s of first order. Therefore consider:
$$\begin{array}{rl}& f:\mathbb{R}\to {\mathbb{R}}^2,\\ & f(t)=\left(\begin{array}{c}\phi (t)\\ \dot{\phi }(t)\end{array}\right).\end{array}$$
Then $\phi$ solves the pendulum equation if and only if $f$ solves
$${\left(\begin{array}{c}u\\ v\end{array}\right)}^{\bullet }=\left(\begin{array}{c}v\\ -\sin (u)\end{array}\right).$$
The energy of the above system is given by $E(t):=\frac{1}{2}{\dot{\phi }}^2–\cos (\phi ).$ $E$ is a constant of motion since:
$$\dot{E}=\dot{\phi }\ddot{\phi }+\dot{\phi }\sin (\phi )=o,$$

The level-sets of the energy of the pendulum. The yellow dots correspond to $E=-1$, the green curves to $-1<E<1$, the blue ones to $E=1$ and the red ones to $E>1$.

if $\phi$ is a solution of the pendulum equation. In the phase space the solutions of the ODE are implicitly given by the level-sets of $E$ :

By the theorem of Picard-Lindelöf there is for every $\left(\begin{array}{c}{u}_0\\ v_0\end{array}\right)\in {\mathbb{R}}^2$ a  unique solution of the initial value problem (IVP) :
$$\begin{array}{rl}& {\left(\begin{array}{c}u\\ v\end{array}\right)}^{\bullet }=X(u,v)=\left(\begin{array}{c}v\\ -\sin (u)\end{array}\right)\\ & \left(\begin{array}{c}u\\ v\end{array}\right)(0)=\left(\begin{array}{c}{u}_0\\ v_0\end{array}\right).\end{array}$$
These solutions lie on the level-sets of $E$.
Now we want to approximate the solution of the IVP by a sequence
$$\dots ,\left(\begin{array}{c}{u}_{-1}\\ v_{-1}\end{array}\right),\left(\begin{array}{c}{u}_0\\ v_0\end{array}\right),\left(\begin{array}{c}{u}_1\\ v_1\end{array}\right),\dots \in {\mathbb{R}}^2$$
such that for the real solution $\left(\begin{array}{c}u(t)\\ v(t)\end{array}\right)$ we have $\underset{p_k}{\underbrace{\left(\begin{array}{c}{u}_k\\ v_k\end{array}\right)}}\approx \left(\begin{array}{c}u\\ v\end{array}\right)(k\delta ),$
where $\delta$ is the time step.

There are several methods to get such a sequence:

1. Forward Euler method:

Let $X:{\mathbb{R}}^2\to {\mathbb{R}}^2$ be the vectorfiled representing the ODE , i.e. $${\left(\begin{array}{c}u\\ v\end{array}\right)}^{\bullet }(k\delta )=X(u(k\delta ),v(k\delta ))=X(p_k).$$

Then with the forward Euler method we get from $p_k$ to $p_{k+1}$ by walking for the time $\delta$ in the direction of $X(p_k)$.

$$p_{k+1}=p_k+\delta X(p_k).$$

The big advantage of this method is, that it is very fast. The disadvantage is illustrated in the following example: For the harmonic oscillator the ODE is given by the vectorfield:
$$X(u,v)=\left(\begin{array}{c}v\\ -u\end{array}\right).$$
Its analytic solution is:
$$\left(\begin{array}{c}u\\ v\end{array}\right)(t)=\left(\begin{array}{cc}\cos (t)& \sin (t)\\ -\sin (t)& \cos (t)\end{array}\right)\left(\begin{array}{c}{u}_o\\ v_0\end{array}\right)$$

The vector field and the corresponding integralcurves of the harmonic oscillator.

The energy is a again a constant of motion: $E(t):=\frac{1}{2}(u^2+v^2)$ and the level-sets of $E$ are circles. If we apply the forward Euler method the solution is a spiral and not a circle. Therefore, the energy gets bigger and bigger and the solution blows up.

The real integralcurve (blue) and the one obtain with the Euler forward method for $\delta =0.5$  (red), together with the corresponding vectorfield (grey).

2. Backward Euler method:

The backward Euler method works in the following way: Find $p_{k+1}$ such that one Euler forward step with respect to the vectorfield $-X$ gives you $p_k$. i.e
$$p_k=p_{k+1}-\delta X(p_{k+1}).$$
That means we have to solve the following non linear system of equations to get from $p_k$ to $p_{k+1}$:
$$p_{k+1}=p_k+\delta X(p_{k+1}).$$

This can for example be done with the classical Newton method. The advantage of the backward Euler method is that the energy does not increase and therefore the solution does not blow up. The disadvantage is that solving a non linear system of equations can be hard and expensive.

The idea behind both Euler methods, that guaranties the convergence to the real solution for $\delta \to 0$ is given by the mean value theorem: For $f\in C^1(\mathbb{R},{\mathbb{R}}^n)$, $p\in \mathbb{R}$ and $\delta >0$ there exist $\tau \in (p,p+\delta )$ such that”
$$f^{\prime }(\tau )=\frac{f(p)-f(p+\delta )}{\delta }.$$

3. Verlet method:

The Verlet method works for equations of the form $$\ddot{f}(t)=F(f(t)),f\in C^2(\mathbb{R},{\mathbb{R}}^n).$$

In order to approximate $f(k\delta )$ by $p_k$ we use the mean value theorem twice:

$$\begin{array}{rl}\ddot{f}(t)& \approx \frac{\dot{f}(t+\frac{\delta }{2})-\dot{f}(t-\frac{\delta }{2})}{\delta }\approx \frac{\frac{f(t+\delta )-f\left(t\right)}{\delta }-\frac{f\left(t\right)-f(t-\delta )}{\delta }}{\delta }\\ & =\frac{f(t+\delta )-2f(t)+f(t-\delta )}{{\delta }^2}.\end{array}$$

This gives us the following recursion formula for the $p_k$’s:
$$\begin{array}{rl}F(p_k)& =\frac{p_{k-1}-2p_k+p_{k+1}}{{\delta }^2}\\ \iff p_{k+1}& =2p_k-p_{k-1}+{\delta }^{2}F(p_k).\end{array}$$

To solve this recursion formula, we start with the initial data $$\left(\begin{array}{c}{p}_0\\ p_1\end{array}\right):=\left(\begin{array}{c}{u}_0\\ v_1\end{array}\right):=\left(\begin{array}{c}f(0)\\ f(\delta )\end{array}\right)$$
with the approximation:
$$\begin{array}{r}\frac{p_0+p_1}{2}\approx f\left(\frac{\delta }{2}\right),\\ \frac{p_1-p_0}{\delta }\approx \dot{f}\left(\frac{\delta }{2}\right),\end{array}$$

this corresponds to the usual initial data of an IVP for a  for a second order ODE. Afterwords we use:

$$\left(\begin{array}{c}{p}_k\\ p_{k+1}\end{array}\right)=\left(\begin{array}{c}{u}_k\\ v_k\end{array}\right)=\left(\begin{array}{c}{v}_{k-1}\\ 2v_{k-1}-u_{k-1}+{\delta }^{2}F(v_{k-1})\end{array}\right).$$

If we consider that as a two dimensional system with $q_k:=\left(\begin{array}{c}{u}_k\\ v_k\end{array}\right)$, we have:
$$q_{k+1}=X(q_k),$$
with:
$$\begin{array}{rl}& X:{\mathbb{R}}^2\to {\mathbb{R}}^2,\\ & X(u,v)=\left(\begin{array}{c}v\\ 2v-u+{\delta }^{2}F(u)\end{array}\right).\end{array}$$

This is a very good choice of discretization, because the  vectorfield $X$ is area preserving:
$$X^{\prime }(u,v)=\left(\begin{array}{cc}0& 1\\ -1& 2+{\delta }^2{F}^{\prime }(u)\end{array}\right)⇝det{X}^{\prime }(u,v)=1.$$
This gives us for $U\subset {\mathbb{R}}^2:$
$$\text{area}(X(U))={\int }_{X(U)}1={\int }_U|det{X}^{\prime }(u,v)|={\int }_{U}1=\text{area}(U).$$

There is a huge theory about area preserving maps and their higher dimensional analog symplectic maps. This maps lead to so called symplectic integrators, that do an excellent job solving ODE’s.

Lets get back to the pendulum equation: $\ddot{\phi }=-\sin (\phi ).$ If we use a clever method for the discretization of $F(\phi )=\sin (\phi )$, the Verlet method gives us a good solution of the pendulum equation. Where the clever discretization of $F$ is the following: $$F(\phi )=\frac{1}{k}\arg (1+k{e}^{i\phi }).$$

Let $M=(V,E,F)$ be an oriented triangulated surface without boundary and $p:V\to {\mathbb{R}}^3$ a realization. In earlier lectures we considered the space of piecewise linear functions on $M$:
$$W_{PL}:=\{\tilde{f}:M\to \mathbb{R}|{\tilde{f}|}_{\sigma }\text{ is affine for all }\sigma \in F\}.$$
A basis of $W_{PL}$ is given by : ${\varphi }_i\in W_{PL}$ defined by ${\varphi }_i(p_j)={\delta }_{ij},$ an  interpolation map is defined by:
$$\begin{array}{rl}{\text{int}}_{PL}:{\mathbb{R}}^V& \to W_{PL},\\ f& \mapsto \tilde{f}=\sum _{i\in V}{f}_i{\varphi }_i.\end{array}$$

The graph of ${\varphi }_i$

For some purposes piecewise constant function are more convenient then piecewise linear ones, therefore, we introduce a new function space on $M$ that consists of piecewise constant functions. In order to do that we assign a domain to each vertex:
For every edge $e_{ij}\in E$ the center is given by $S_{ij}=\frac{1}{2}(p_i+p_j)$ analog  for every triangle $\sigma =\{i,j,k\}\in F$ the center of mass is defined as $S_{ijk}=\frac{1}{3}(p_i+p_j+p_k)$.

If we connect the center of the edges with the center of mass the triangle is divided into three parts with equal area. Now we can assign to each vertex $p_i$ the parts of the adjacent triangles that contain $p_i$ and denote it with $D_i$.

The domain $D_i$ (gray colored) associated to the vertex $p_i$

For the area of $D_i$ we obtain:
$$A_i:=\frac{1}{3}\sum _{\sigma =\{i,j,k\}\in F}{A}_{\sigma }.$$

The space of piecewise constant functions is defined in the following way:
$$W_{PC}:=\{\hat{f}:M\to \mathbb{R}|{\hat{f}|}_{D_i}\text{ is constant for all }i\in V\}.$$
Also on $W_{PC}$ we have a canonical basis ${\mathrm{\Psi }}_i\in W_{PC}$ defined by ${\mathrm{\Psi }}_i(p_j)={\delta }_{ij}$ and an interpolation map:
$$\begin{array}{rl}{\text{int}}_{PC}:{\mathbb{R}}^V& \to W_{PC},\\ f& \mapsto \hat{f}=\sum _{i\in V}{f}_i{\mathrm{\Psi }}_i.\end{array}$$

Graph of ${\mathrm{\Psi }}_i$

There is a similar construction that assigns to each vertex an face and works purely combinatorial: For a triangulated $M$ surface one can define its discrete dual surface $M^{\ast }=(V^{\ast },E^{\ast },F^{\ast })$. The faces of $M$ correspond to the vertices of $M^{\ast }$ and the vertices of $M$ to the faces of $M^{\ast }$. The set of edges stays the same. Two vertices in $M^{\ast }$ are connected by  an edge if the corresponding faces in $M$ are adjacent. For more informations see here.

A piece of an triangulated surface $M$ and its dual $M^{\ast }$. The red vertices and edges belong to $M$ and the blues ones to $M^{\ast }$.

Since in general it is not easy to compute the area of the dual faces with the metric of the original one, we used the above construction of the $D_i$’s to define $W_{PC}$. On $W_{PL}$ we defined the usual scalar product $⟨⟨\tilde{f},\tilde{g}⟩{⟩}_{PL}:={\int }_M\tilde{f}\tilde{g}$ and computed the matrix $B$ representing it with respect to the basis $\{{\varphi }_1,\dots {\varphi }_N\},$ i.e. $b_{ij}:=⟨⟨{\varphi }_i,{\varphi }_j⟩{⟩}_{PL}.$ This gave us a scalar product on ${\mathbb{R}}^V$:

$$⟨⟨f,g⟩{⟩}_{PL}=f^{t}Bg=\sum _{i,j\in V}{f}_i{g}_j{b}_{ij}.$$

On the same way we define an scalarproduct on $W_{PC}$:
$$\begin{array}{rl}⟨⟨\cdot ,\cdot ⟩{⟩}_{PC}& :W_{PC}\times W_{PC}\to \mathbb{R},\\ ⟨⟨\hat{f},\hat{g}⟩{⟩}_{PC}& :={\int }_M\hat{f}\hat{g}.\end{array}$$
The matrix $A$ representing it with respect to the basis  $\{{\mathrm{\Psi }}_1,\dots {\mathrm{\Psi }}_N\}$ is a diagonal matrix since:
$$a_{ij}:=⟨⟨{\mathrm{\Psi }}_i,{\mathrm{\Psi }}_j⟩{⟩}_{PC}={\int }_M{\mathrm{\Psi }}_i{\mathrm{\Psi }}_j=A_i{\delta }_{ij}.$$
This gives us a new  scalar product on ${\mathbb{R}}^V:$
$$⟨⟨f,g⟩{⟩}_{PC}:=f^{t}Ag=\sum _{i\in V}{f}_i{g}_i{a}_{ii}.$$
Note that the basis $\{{\mathrm{\Psi }}_1,\dots {\mathrm{\Psi }}_N\}$ is orthogonal with respect to $⟨⟨\cdot ,\cdot ⟩{⟩}_{PC}$, while the same is not true for $\{{\varphi }_1,\dots {\varphi }_N\}$ with respect to $⟨⟨\cdot ,\cdot ⟩{⟩}_{PL}$.
Now we have two scalar products on ${\mathbb{R}}^V$ one represented by a symmetric positive definite matrix $B$ and one by a positive definite diagonal one $A$.
Our aim is to define a discrete laplace operator $\mathrm{\Delta }:{\mathbb{R}}^V\to {\mathbb{R}}^V$ analog to the smooth laplacian we discussed earlier. We will see that the discrtelaplacian depends on the choice of scalarproduct.
On a smooth surface $M$ we have:
$$⟨⟨\mathrm{\Delta }f,g⟩⟩={\int }_M\mathrm{\Delta }fg=–{\int }_M⟨\text{grad}f,\text{grad}g⟩.$$
This identity should also be true in the discrete case and will serve us as definition for the discrete laplacian. We already defined the bilinear operator
$$\begin{array}{rl}L:{\mathbb{R}}^V\times {\mathbb{R}}^V& \to \mathbb{R},\\ f,g& \mapsto –{\int }_M⟨\text{grad}f,\text{grad}g⟩,\end{array}$$

and discussed the relation between  quadratic forms, symmetric bilinear maps and maps between dual spaces ( lecture: dirichlet energy 2 ).
$$\begin{array}{lcl}{W}_{PL}& \stackrel{E_D}{⟶}& W_{PL}^{\ast }\\ \uparrow \text{interpol.}& & \downarrow {\text{interpol.}}^{\ast }\\ {\mathbb{R}}^V& \stackrel{L}{⟶}& {\mathbb{R}}^{V\ast }\end{array}$$
In the finite case the quadratic form, bilinear form and map between the dual spaces can all be represented by the same matrix $L$ .

For $f,g\in {\mathbb{R}}^V$ we have:
$$f^{t}Lg=-{\int }_M⟨\text{grad}f,\text{grad}g⟩\stackrel{!}{=}{\int }_M\mathrm{\Delta }fg=⟨⟨\mathrm{\Delta }f,g⟩⟩.$$
This gives us one laplace operator for each scalarproduct on ${\mathbb{R}}^V$:
$$\begin{array}{r}{f}^{t}Lg=\{\begin{array}{c}{\left({\mathrm{\Delta }}_{PL}f\right)}^{t}Bg,\\ {\left({\mathrm{\Delta }}_{PC}f\right)}^{t}Ag,\end{array}\\ \Rightarrow \{\begin{array}{c}{\mathrm{\Delta }}_{PL}={\left(L{B}^{-1}\right)}^t=B^{-1}L,\\ {\mathrm{\Delta }}_{PC}={\left(L{A}^{-1}\right)}^t=A^{-1}L.\end{array}\end{array}$$

In general it is complicated and expensive to compute $B^{-1}$, such that it is much easier to use ${\mathrm{\Delta }}_{PC}$ instead of ${\mathrm{\Delta }}_{PL}$. Witch of the laplacians one has to choose depends on the nature of the problem.

One famous problem form physics and geometry is the Poisson problem:
Given $\rho \in {\mathbb{R}}^V$ find $u\in {\mathbb{R}}^V$ such that:
$$\mathrm{\Delta }u=-\rho .$$
For Poisson problems both laplacians can be used since:

$${\mathrm{\Delta }}_{PL}u={\left(L{B}^{-1}\right)}^{t}u=-\rho \iff Lu=-B\rho .$$

Example: Suppose $M$ consists of an uniform conducting material and $\rho :M\to \mathbb{R}$ is a charge density. For the electric field $E$ on $M$ there holds:
$$\text{div}E=\rho ,$$
and for the potential $u:M\to \mathbb{R}$ we have
$$E=-\text{grad}u.$$Taking this together we have to solve:
$$\mathrm{\Delta }u=-\rho ,$$
in order to determine the potential for a given density.
In this setting piecewise constant functions are much more convenient to describe the charge density then piecewise linear ones. Therefore we use ${\mathrm{\Delta }}_{PL}$ to solve this problem.

Let $\mathrm{\Sigma }=(V,E,F)$ be a triangulated surface (with boundary). A realization of the surface in ${\mathbb{R}}^3$ is given by a map $p:V\to {\mathbb{R}}^3$ such that $p_i,p_j,p_k$ form a non degenerated triangle in ${\mathbb{R}}^3$ for all $\{i,j,k\}\in \mathrm{\Sigma }$, (i.e. $p_i,p_j,p_k$ do not lie on a line). Some problems, like the plateau problem do not depend on the position of the vertices in the space, but only on the length of the edges:
$$l_{ij}=l_{ji}=|e_{ji}|i,j\in E.$$
If we have a realization of the surface in ${\mathbb{R}}^3$ there holds: $l_{ij}=l_{ji}=|p_i-p_j|$ and the edge length’s satisfy the triangle inequalities:
$$\text{For }\{i,j,k\}\in F\Rightarrow \{\begin{array}{c}{l}_{ij}<l_{jk}+l_{ki},\\ l_{ki}<l_{ij}+l_{jk},\\ l_{jk}<l_{ki}+l_{ij}.\end{array}$$

Definition: A metric on a triangulated surface assigns  a length to each oriented edge such that the triangle inequalities are satisfied.

Note that every realization $p:V\to {\mathbb{R}}^3$ induces a metric on the underlying triangulated surface. On the other hand a triangulated surface with metric has several realization in ${\mathbb{R}}^3$, such that the metrics coincide. On a triangulated surface with metric every triangle has the structure of an euclidean triangle in ${\mathbb{R}}^2$, that allows us to interpolate $f:V\to \mathbb{R}$ to a piecewise affine function $\tilde{f}:C\left(\mathrm{\Sigma }\right)\to \mathbb{R}$ and define the Dirichlet energy $E_D(f)$ as in the last lectures.

 In particular, we can define angles and area using the cosine theorem:
$$\cos ({\alpha }_i)=\frac{l_{jk}^2-l_{ji}^2-l_{ki}^2}{2l_{ki}{l}_{ij}}$$
By the triangle inequalities we obtain that $0<\gamma <\pi$ and the angles are well defined. With $\sin ({\alpha }_i)=\sqrt{1-\cos ({\alpha }_i{)}^2}$ we can define the cotangents weights and the area of the triangle $\sigma =\{i,j,k\}$ :
$$A_{\sigma }=\frac{1}{2}\sin ({\alpha }_i)l_{ij}{l}_{ki}.$$
For a realization $p$ of $\mathrm{\Sigma }$ with vertex positions $p_i,p_j,p_k$ and edges $e_{ij}=p_j-p_i,\dots$ we obtain for $\sigma =\{i,j,k\}\in F$:
$$A_{\sigma }(p)=\frac{1}{2}|(p_j-p_i)\times (p_k-p_i)|=\frac{1}{2}|e_{ij}\times e_{ik}|=\frac{1}{2}\sin ({\alpha }_i)\underset{=l_{ij}}{\underbrace{|e_{ij}|}}\underset{=l_{ki}}{\underbrace{|e_{ki}|}},$$
and the total area is give by:
$$A(p)=\sum _{\sigma \in F}{A}_{\sigma }(p).$$

Definition: The following problem is the Plateau problem: Let $\mathrm{\Sigma }$ be a triangulated surface with boundary and $q:V^{\partial }\to {\mathbb{R}}^3$ a prescribed function. We are looking for a realization $p:V\to {\mathbb{R}}^3$ of $\mathrm{\Sigma }$ such that:

1.  ${p|}_{V^{\partial }}=q,$
2.  $p$ has the smallest area among all realizations $\tilde{p}$ of $\mathrm{\Sigma }$, with ${\tilde{p}|}_{V^{\partial }}=q$.

In order to be able to solve the  Plateau problem we have to consider variations of $p$:

Definition: A variation of $p$ with constant boundary assigns to each vertex  $i\in V$ a curve $t\mapsto p_i(t),$ $t\in (-ϵ,ϵ),$ such that:

1. $p_i(0)=p_i$ for all $i\in V.$
2. $p_j(t)=q_j$ for all $t\in (-ϵ,ϵ)$ if $j\in V^{\partial }$.

Definition: $p$ is called a minimal surface if it is a critical point of the area functional, i.e. :
$${\frac{\partial }{\partial t}|}_{t=0}A(p(t))=0,$$
for all variations of $p$ with constant boundary.

Lemma: $p$ is a minimal surface if and only if For each $i\in \stackrel{\circ }{V}$ and each curve $t\mapsto p_i(t),$ $t\in (-ϵ,ϵ)$ we have:
$${\frac{\partial }{\partial t}|}_{t=0}A(p_i(t))={\frac{\partial }{\partial t}|}_{t=0}\sum _{\sigma \in F|i\in \sigma }{A}_{\sigma }(p_i(t))=0.$$

I.e. $p$ is a critical point of the area functional with respect to all variations with constant boundary, if and only if it is a critical point for all those variations that only move one interior vertex.

Proof: $”\Rightarrow ”:$ clear.
$”\Leftarrow :”$ Similar to the following fact from analysis:
All directional derivatives of a function $f:U\to \mathbb{R},U\subset {\mathbb{R}}^n$ vanish at $p\in U$ if and only if all partial derivatives of $f$ vanish at $p$.

Theorem: $p$ is a minimal surface if and only if all tree components of $p$ are discrete harmonic functions with respect to the metric induced on $\mathrm{\Sigma }$ by $p$.

Proof: Fix $i\in \stackrel{\circ }{V}$ and let  $t\mapsto p_i(t),$ $t\in (-ϵ,ϵ)$ be a variation, that only moves $p_i$ and $Y:=p_i^{\prime }(o).$ Further we define $\hat{F}:=\{\{i,j,k\}\in F|\text{ positive oriented}\}$ to be the set of positive oriented triangles that contain the vertex $i$. We have to show that there holds:
$$\begin{array}{rl}0={\frac{\partial }{\partial t}|}_{t=0}A(p_i(t))& ={\frac{\partial }{\partial t}|}_{t=0}\sum _{\{i,j,k\}\in \hat{F}}|(p_j-p_i(t))\times (p_k-p_i(t))|\\ & ={\frac{\partial }{\partial t}|}_{t=0}\sum _{\{i,j,k\}\in \hat{F}}|e_{ij}(t)\times e_{ik}(t)|.\end{array}$$

Before we start with the proof remember these identities from analysis$\mathrm{\setminus }$linear algebra

1.  For a function $t\to v(t)\in {\mathbb{R}}^3$ there holds:
$$\begin{array}{r}|v{|}^{\prime }=\frac{⟨v,v^{\prime }⟩}{|v|}\end{array}$$
2. For $a,b,c,d\in {\mathbb{R}}^3$ we have:
$$\begin{array}{r}⟨a\times b,c\times d⟩=det\left(\begin{array}{cc}⟨a,c⟩& ⟨a,d⟩\\ ⟨b,c⟩& ⟨b,d⟩\end{array}\right)=⟨a,c⟩⟨b,d⟩-⟨a,d⟩⟨b,c⟩.\end{array}$$
With $e_{i,j}=p_j-p_i$ we get $\dot{e_{ij}}=-Y$ and can compute now the variation of the area functional:
$$\begin{array}{rl}{\frac{\partial }{\partial t}|}_{t=0}& A(p_i(t))=\frac{1}{2}\sum _{\{i,j,k\}\in \hat{F}}\frac{⟨e_{ij}\times e_{ik},-Y\times e_{ik}–e_{ij}\times Y⟩}{|e_{ij}\times e_{jk}|}\\ & =\sum _{\{i,j,k\}\in \hat{F}}\frac{-⟨e_{ij},Y⟩|e_{ik}{|}^2+⟨e_{ik},Y⟩⟨e_{ij},e_{ik}⟩-⟨e_{ik},Y⟩|e_{ij}{|}^2+⟨e_{ij},Y⟩⟨e_{ij},e_{ik}⟩}{2|e_{ij}\times e_{jk}|}\\ & =⟨Y,\sum _{\{i,j,k\}\in \hat{F}}\frac{-e_{ij}|e_{ik}{|}^2-e_{ik}|e_{ij}{|}^2+⟨e_{ij},e_{ik}⟩(e_{ij}+e_{ik})}{2|e_{ij}\times e_{jk}|}⟩\\ & =⟨Y,\sum _{\{i,j,k\}\in \hat{F}}\frac{(-⟨e_{ik},e_{ik}⟩+⟨e_{ij},e_{ik}⟩)e_{ij}+(-⟨e_{ij},e_{ij}⟩+⟨e_{ij},e_{ik}⟩)e_{ik}}{2|e_{ij}\times e_{jk}|}⟩.\end{array}$$
With $$-⟨e_{ik},e_{ik}⟩+⟨e_{ij},e_{ik}⟩=⟨e_{ik},e_{ij}-e_{ik}⟩=⟨e_{ik},e_{kj}⟩,$$
and
$$-⟨e_{ij},e_{ij}⟩+⟨e_{ij},e_{ik}⟩=⟨e_{ij},e_{ik}-e_{ij}⟩=⟨e_{ij},e_{jk}⟩,$$
We have:
$${\frac{\partial }{\partial t}|}_{t=0}A(p_i(t))=⟨Y,\sum _{\{i,j,k\}\in \hat{F}}\left(\frac{⟨e_{ik},e_{kj}⟩e_{ij}}{2|e_{ij}\times e_{jk}|}\right)+\sum _{\{i,j,k\}\in \hat{F}}\left(\frac{⟨e_{ij},e_{jk}⟩e_{ik}}{2|e_{ij}\times e_{jk}|}\right)⟩.$$

In an earlier lecture  we showed that: $\cot ({\beta }_{ik})=\frac{⟨-e_{ij},e_{jk}⟩}{|e_{ij}\times e_{jk}|}$ and $\cot ({\alpha }_{ik})=\frac{⟨e_{il},e_{kl}⟩}{|e_{ik}\times e_{il}|}.$

With an index shift in the first sum $(i,j,k)\to (i,k,l)$ we have now : $$\begin{array}{rl}{\frac{\partial }{\partial t}|}_{t=0}A(p_i(t))& =⟨Y,-\sum _{\{i,j,k\}\in \hat{F}}\frac{1}{2}(\cot ({\beta }_{ik})+\cot ({\alpha }_{ik}))e_{ik}⟩\\ & =⟨Y,-\sum _{\{i,k\}\in \hat{E}}\frac{1}{2}(\cot ({\beta }_{ik})+\cot ({\alpha }_{ik}))(p_k–p_i)⟩.\end{array}$$

Since this has to hold for all variations $Y$ of $p_i$ and for all $i\in \stackrel{\circ }{V}$we finally get:
$${\frac{\partial }{\partial t}|}_{t=0}A(p(t))=0\iff \sum _{\{i,k\}\in \hat{E}}\frac{1}{2}(\cot ({\beta }_{ik})+\cot ({\alpha }_{ik}))(p_k–p_i)=0\mathrm{\forall }i\in \stackrel{\circ }{V}.$$

This means $p$ is an minimal surface if and only if all components of $p$ are harmonic with respect to the metric on $\mathrm{\Sigma }$ induced by $p$.

Let $M=(V,E,F)$  be a triangulated domain in the plane where $V$ denotes the set of vertices, $E$ the set of edges and $F$ the set of triangles.  We consider the set of functions on the vertices  :
$$\begin{array}{r}{\mathbb{R}}^V:=\left\{\begin{array}{cc}f:& V\to \mathbb{R}\\ & i\mapsto f_i\end{array}\right\}.\end{array}$$
By interpolation to a piecewise linear function we get:  $\hat{f}:M\to \mathbb{R}$.
$${\mathbb{R}}^V\stackrel{\text{interpol.}}{\to }{W}_{PL}:=\{\hat{f}:M\to \mathbb{R}|{\hat{f}|}_{\sigma }\text{ is affine for all }\sigma \in F\}.$$
In the last lecture the Dirichlet energy of  $\hat{f}\in W_{PL}$ was defined as:
$$E_D(\hat{f})=\frac{1}{2}{\int }_M{\left|\text{grad}\tilde{f}\right|}^2=\frac{1}{4}\sum _{(i,j)\in \hat{E}}(\cot ({\alpha }_{ij})+\cot {\beta }_{ij}){(f_i-f_j)}^2.$$
$E_D:W_{PL}\to \mathbb{R}$ is a quadratic form, which can be represented by matrix $A\in {\mathbb{R}}^{V\times V}:$
$$E_D(f)=(f_1,\dots ,f_N)\left(\begin{array}{ccc}{a}_{11}& \dots & a_{1N}\\ ⋮& \ddots & ⋮\\ a_{N1}& \dots & a_{NN}\end{array}\right)\left(\begin{array}{c}{f}_1\\ ⋮\\ f_N\end{array}\right),$$
where $N:=|V|$ and
$$a_{ij}:=\{\begin{array}{rcl}-\frac{1}{4}(\cot ({\alpha }_{ij})+\cot {\beta }_{ij})& \text{for}& i\ne j,\{i,j\}\in E,\\ \frac{1}{4}\sum _{(i,k)\in \hat{E}}(\cot ({\alpha }_{ik})+\cot {\beta }_{ik})& \text{for}& i=j,\\ 0& \text{else}.& \end{array}$$
One can easily check that: $A\left(\begin{array}{c}1\\ ⋮\\ 1\end{array}\right)=0$ and even further there holds: $\text{ker }A=\mathbb{R}\left(\begin{array}{c}1\\ ⋮\\ 1\end{array}\right).$