Trying out the Schrödinger

Time-independent Schrödinger equation constitutes the following linear eigenvalue problem: find $u : \mathbb{R}^3 \rightarrow \mathbb{C}$ and $\lambda \in \mathbb{C}$ such that
$$-\frac{h^2}{4\pi \mu} \Delta u + V(\boldsymbol{x}) u = \lambda u,$$ where $h$ and $\mu$ are some constants and $V(\boldsymbol{x})$ is a given potential. For hydrogen atom we have
$$V(r) = -\frac{e^2}{4\pi \epsilon_0 r},$$ where $e$ and $\epsilon_0$ are some additional constants and $r$ is the distance from the nucleus. Writing all the constants as one constant
$$C = - \frac{\mu e^2}{\epsilon_0 h^2}$$ and rescaling the eigenvalue leads to
$$-\Delta u - C r^{-1} u = \widetilde{\lambda} u.$$ Let $\Omega$ be the unit ball $B(0,1)$ and consider the homogeneous Dirichlet boundary condition $u|_{\partial \Omega} = 0$. Weak formulation for this problem is: find $u \in H^1_0(\Omega)$ such that
$$(\nabla u, \nabla v)-C(r^{-1}u,v)=\widetilde{\lambda} (u,v)$$ holds for every $v \in H^1_0(\Omega)$.

Setting $C=1$ and solving this with piecewise linear elements on my Matlab FEM code gives eigenvectors whose midsections look as follows:

First 12 lowest eigenvectors of the Schrödinger equation with constants normalized to 1.

You may compare them to the solutions given at Wikipedia. Can you find them all?

On ray marching

One way to define a surface in three-dimensions is through a function $f:\mathbb{R}^3\rightarrow\mathbb{R}$ which satisfies
$$f(\boldsymbol{x})=0,$$ for every point $\boldsymbol{x}\in \mathbb{R}^3$ that belongs to the surface. Thus, the surface $S$ can be defined as
$$S=\{\boldsymbol{x} \in \mathbb{R}^3\,|\,f(\boldsymbol{x})=0\}.$$ For one to visualize the surface $S$, the intersection point $\boldsymbol{p} \in \mathbb{R}^3$ of a ray---starting from the camera position $\boldsymbol{c} \in \mathbb{R}^3$ to some direction $\boldsymbol{d} \in \mathbb{R}^3$, $\|\boldsymbol{d}\|=1$---and the surface must be sought. After these intersection points are known for sufficiently many rays, a visualization can be constructed (some examples follow).

One way to find the intersection point $\boldsymbol{p}$ is through analytical calculations: the surface defined by $S$ is sometimes so simple that this calculation can be done by hand. An example that considers the intersection of a ray and a sphere can be found from Wikipedia. This is exactly what the traditional ray tracing is all about.

In a more general setting, no analytical relationship exists for finding the intersection point $\boldsymbol{p}$. In such a case, one must rely to numerical methods.

A proficient applied mathematician would immediately attempt to apply, e.g., Secant method to the problem: find $t \in \mathbb{R}$ such that
$$f(\boldsymbol{c}+t\boldsymbol{d}) \approx 0.$$ There is, however, a caveat. For a general $f$ the function
$$g(t):=f(\boldsymbol{c}+t\boldsymbol{d})$$ might have several roots and in Secant method (or similar) nothing guarantees that a root nearest to the camera is recovered. You might end up drawing the wrong side of the object of interest! Thus, the applied numerical method should be guaranteed to find a root nearest to the camera.

A somewhat trivial method to find the nearest (or approximately nearest) root is to fix a sufficiently small constant step size $\Delta t$ and then find the smallest $n \in \mathbb{N} \cup \{0\}$ such that
$$f(\boldsymbol{c}+n \Delta t \boldsymbol{d})\approx 0.$$ This is known as ray marching and is indeed a valid method for visualizing the surface $S$. The problem is that it's tremendously slow since one must try every possible $n$ starting from zero!

Fortunately, there exists a clever trick to make the computation faster by ensuring that $f$ defines a distance field of $S$. This means that at every point $\boldsymbol{x}$ in the space $\mathbb{R}^3$ the value $f(\boldsymbol{x})$ represents the minimum distance from $\boldsymbol{x}$ to the surface $S$. This sounds like a hard restriction to satisfy but for many fractals the distance fields (or approximate distance fields) have already been derived by enthusiasts.

The algorithm goes as follows:

 1. $\boldsymbol{p}:=\boldsymbol{c}$
 2. $\boldsymbol{p}:=\boldsymbol{p}+f(\boldsymbol{p})\boldsymbol{d}$
 3. If $f(\boldsymbol{p})$ is larger than some tolerance go back to step 2. Otherwise, break.

Output of the prescribed algorithm is the intersection point $\boldsymbol{p}$.

If you happen to own a sufficiently fast graphics processing unit and know the function $f$ for some fractal then this method can be implemented in real time using, e.g., OpenGL Shading Language. See the following Shadertoy entry I constructed as an example:

Visualizing a mandelbulb fractal. See https://www.shadertoy.com/view/MtX3zr.

Getting the coloring and other visuals right requires of course lots of experience (or trial-and-error) but this is the basic idea for generating such animations. I tried to add some comments to the shader code and attempted to use same symbols as in this post.

P.G. Ciarlet on the words "linear triangle" and "biquadratic rectangle"

Since the readers of this article are perfectly aware that a triangle may be isoceles, but certainly not genuinely "linear", or that a rectangle may be square, but certainly not genuinely "biquadratic", I have blithely committed these serious abus de langage, which clearly convey more information than the names that I had originally chosen.

P.G. Ciarlet, Handbook of Numerical Analysis, Part 1

Exporting to vectorized PDF images in MATLAB

Keywords: painters, pdf, trisurf, white lines, patch

Ever tried to export nice looking, smoothly shaded pdf figures from MATLAB and ended up with something like this?

Would you prefer something like this?

If the answer to both questions is yes, then read on.

A plenty of forum posts in Mathworks' website propose different kinds of magic tricks to solve this problem. However, none of the tricks seemed to work for my case where I was drawing a surface with trisurf and exporting it to vectorized PDF with the following code.

trisurf(mesh.t',mesh.p(1,:)',mesh.p(2,:)',...
                  ufun(mesh.p(1,:),mesh.p(1,:))');
shading interp;
print -painters -dpdf -r600 vel.pdf

This always ended up with the annoying white lines as demonstrated in the first figure. Finally I realized that trisurf uses underneath the patch command to actually draw the triangles. I came up with a simple fix of setting up interpolated edge colors to every triangle drawn by the patch command as follows:

H = trisurf(mesh.t',mesh.p(1,:)',mesh.p(2,:)',...
                  ufun(mesh.p(1,:),mesh.p(1,:))');
set(H,'EdgeColor','interp');
shading interp;
print -painters -dpdf -r600 vel.pdf

This gives the beatiful plot in the second figure. For cropping the whitespace around the pdf I propose the usage of the linux command line tool pdfcrop.

3D progress

Today I added preliminary support for 3D bilinear and linear forms using tetrahedral P1 elements. Unfortunately, the visualization is based on MATLAB using MATLAB.jl Julia package (https://github.com/lindahua/MATLAB.jl) and the mesh generation is also based on a MATLAB script package iso2mesh (http://iso2mesh.sourceforge.net/cgi-bin/index.cgi) which is again based on tetgen (http://wias-berlin.de/software/tetgen/).

Thus, I haven't yet figured out how I should package the code properly for distribution. Meanwhile, here is a screenshot of solving the Poisson equation
$$\Delta u = 1$$ on a unit cube $[0,1]^3$ with zero Dirichlet boundary conditions and the corresponding mesh:

Cross section at $z=0.5$.

The mesh. Color indicates simply the $z$-coordinate.

Improved assembly

Yesterday I learned how to do variable argument macro's in Julia and based on my new knowledge implemented the bilinear assembly (and linear assembly) with material parameters using them. Thus, now there is only a single macro for bilinear assembly and bilinear assembly with a material parameter (or some other node-wise defined parameter). Next goals for the Julia-code:

• Assembly of the boundary terms
• Faster mesh refinement
• Simple geometry representations for 2D-domains
• 3D-linear tetrahedral element
• ...
• and in a distant future
• ...
• Geometry tools
• Additional elements: p-refinement, vectorial elements (i.e. Raviart-Thomas)
• Isoparametric mappings
• Simple geometry representations for 3D-domains

Some of these might be too much to realize in practice but there is always hope.

Progress on Julia code

Today I implemented a not-so-well-thought-out way to have node-wise defined material parameters in the bilinear forms. The whole assembly process is implemented as a macro, since that way there won't be any performance losses related to e.g. evaluation of some anonymous function defining the bilinear form. Such approach however forces me to build an additional macros for different bilinear forms with different amount of nodal parameters. I haven't yet come up with any particularly elegant solution and because of that I currently just have two different macros like

macro bilin_asm(bilin_form,mesh)
...
end

macro bilin_asm_mat(bilin_form,mesh,mat)
...
end

which are then used as follows:

S = @bilin_asm(ux*vx+uy*vy,mesh)

S = @bilin_asm_mat(a*(ux*vx+uy*vy),mesh,mat)

Here mat is a vector of size $N \times 1$ where $N$ is the amount of nodes in the mesh. Using this approach I can quite easily solve some simple non-linear Poisson equations using e.g. Picard iteration. As an example I decided to try the minimal surface problem: Find the function $u$ for which $$\nabla \cdot \left(\frac{\nabla u}{\sqrt{1+|\nabla u|^2}}\right) = 0$$ on the domain $\Omega = (0,1)^2$ and $u = g$ on the boundary $\partial \Omega$. Such function has minimal surface area of all functions satisfying the boundary condition $u|_{\partial \Omega} = g$. The weak formulation reads: Find $u \in V := \{ w \in H^1(\Omega)\,|\,w|_{\partial \Omega} = g\}$ such that $$\int_\Omega \frac{\nabla u \cdot \nabla v}{\sqrt{1+|\nabla u|^2}} \, \mathrm{d}x = 0$$ for every $v \in H^1_0(\Omega)$. The linearized problem using the Picard iteration reads: Find $u_k \in V$ such that $$\int_\Omega \frac{\nabla u_k \cdot \nabla v}{\sqrt{1+|\nabla u_{k-1}|^2}} \, \mathrm{d}x = 0$$ for every $v \in H^1_0(\Omega)$. Initial guess $u_0$ can be for example some constant. In my test case I have $$g = \sin(2 \pi x - 7/10) \sin(2\pi y-3/10)$$ and the initial guess $u_0(x) = g(x)$.

Initial guess

 

First iteration

Second iteration

Julia for finite elements, part 2

A small update to the previous post. I'm starting to approach the vectorized MATLAB speed. The new results are as follows.

(The highest DOF case was omitted in Julia because my crappy 'n slow 5-minute mesh refiner hogged all the memory.) Some tricks to speed up the code:

• Devectorize, devectorize, devectorize, ... Write everything that you possibly can using loops. Even some simple matrix-vector-multiplications if you do them inside loops.
• Do not use too much array slicing or anonymous functions (or reduntant function calls) inside loops. This will slow things down since function calls seem to have some overhead.
• Preallocate arrays if possible. Otherwise use 1D arrays and push!-command.

I'm happy with the results of my test: Expect some Julia-sorcery in the future! The next target is parallelization.

Julia for finite elements

There has been lots of hype in the internet regarding Julia, JIT-compiled programming language for scientific computing. The surprising thing is that the compiler doesn't expect that the code is vectorized and quite often the devectorized versions are even faster. This weekend I ran some preliminary tests of 2D finite element bilinear form assembly code.

My Julia knowledge is very limited and I have learned everything during this weekend. However, still I get very promising improvements over the for-looped MATLAB code that we use for teaching finite elements at my university. The vectorized MATLAB code has been under development for 10 years and a similar code is used for research. With some tweaking of Julia code I'd probably be able to close the gap between the blue and red lines. As a conclusion, Julia project is definitely something that scientific computation people should be following closely.

Poisson equation by the finite element method

(If the reader is interested in the details, she/he is advised to read an excellent and cheap book Numerical Solution of Partial Differential Equations by the Finite Element Method by Claes Johnson as a reference.)

Let us consider the boundary value problem: Find $u \in H^1_0(\Omega)$ such that
$$
\begin{aligned}
-\Delta u &= F,\quad\text{in $\Omega$},\\
u &= 0,\quad\text{on $\partial \Omega$},\end{aligned}$$ where $\Omega = (0,1)^2$. Analytical solution exists for some $F$ but if $\Omega$ is arbitrary, the typical approach is to use the finite element method: The problem is casted into the weak formulation and then solved by choosing an appropriate discrete subspace $V_h \subset H^1_0(\Omega)$ and by testing with all the basis functions of $V_h$ to get a matrix equation. The weak formulation reads: Find $u \in V:=H^1_0(\Omega)$ such that
$$(\nabla u, \nabla v) = (F,v),$$for every $v \in V$ where we use the notation$$(a,b)_\Omega := \int_\Omega a\cdot b \, \mathrm{d}x = (a,b).$$After we have defined a basis $\{\phi_j\}, j=1,\dots,N$ for the discrete space $V_h \subset V$ (of course $V_h = \text{span}(\phi_1, \dots \phi_N)$) we get the discrete weak formulation: Find $u_h \in V_h$ such that
\begin{equation}
\label{eq:discreteweak}
(\nabla u_h, \nabla v_h) = (F,v_h),
\end{equation} for every $v_h \in V_h$. Céa's lemma guarantees that the solution $u_h$ is the best in the space $V_h$ as long as $V_h \subset V$. Furthermore, if we define the basis $\phi_j, j=1,\dots,N$ by using a triangulation $\mathcal{T}_h = \{K\}$ ($h$ refers to the maximum circumcircle diameter of $K$'s) of the domain $\Omega = \cup_{K \in \mathcal{T}_h} \overline{K}$ such that
$$\phi_j(N_i) = \delta_{ij}$$ where $\delta_{ij}$ is the Kronecker delta symbol and such that $\phi_j \in \mathcal{P}_1(K)$ for every $j$ then by using something known as the Bramble-Hilbert lemma and some scaling argument we can prove that
$$\|u-u_h\|_{H^1(\Omega)} \leq C h |u|_{H^2(\Omega)},$$ if the domain $\Omega$ contains no "evil singularities" and the load $F$ is smooth enough. Such a priori estimate gives some sort of guarantee that the numerical solution $u_h$ will eventually start to be twice as close to the real solution if we refine the mesh such that $h$ gets halved. This is something we next verify numerically by building our own solver.

By substituting $u_h = \sum_{j=1}^N a_j \phi_j$ to the equation \eqref{eq:discreteweak} and by testing with just the basis functions we get
$$ \sum_{j=1}^N a_j (\nabla \phi_j, \nabla \phi_i) = (F,\phi_i),\quad i=1,\dots,N$$which can be written by using matrix notation as
$$ Sx = f,$$where $S_{ij} = (\nabla \phi_i, \nabla \phi_j)$, $f_i = (F,\phi_i)$ and $x_i = a_i$. Calculation of the entries $S_{ij}$ is rather cumbersome in the global coordinate frame and therefore we will do the numerical integration on a reference triangle $\hat{K}$ with vertices $(0,0), (1,0), (0,1)$.

Exercise: Suppose we have a global triangle $K$ with vertices $n_1,n_2$ and $n_3$. Then the affine mapping $\xi: \hat{K} \rightarrow K$ can be written explicitly as
\begin{equation}
\label{eq:affmap}
\xi(\hat{x}) = A\hat{x}+b,
\end{equation} where
$$A = \begin{bmatrix}n_2 -n_1 &n_3-n_1\end{bmatrix},\quad b=n_1,$$ and $n_j$ for all $j$ are considered as column vectors. $\blacksquare$

To make the calculation more efficient, we note that in the summation
\begin{equation}
\label{eq:summand}
S_{ij} =  (\nabla \phi_i, \nabla \phi_j)_\Omega = \sum_{K \in \mathcal{T}_h} (\nabla \phi_i, \nabla \phi_j)_K,
\end{equation} most of the summands integrate to zero. Thus, we will end up with a more efficient procedure if we're not even trying to evaluate such terms in the first place. It is quite straightforward to see that if we just calculate the matrices
\begin{equation} \label{eq:localstiffness} \hat{S}_K = \begin{bmatrix}
(\nabla \phi_{K_1},\nabla \phi_{K_1})_K & (\nabla \phi_{K_1},\nabla \phi_{K_2})_K & (\nabla \phi_{K_1},\nabla \phi_{K_3})_K \\
(\nabla \phi_{K_2},\nabla \phi_{K_1})_K & (\nabla \phi_{K_2},\nabla \phi_{K_2})_K & (\nabla \phi_{K_2},\nabla \phi_{K_3})_K \\
(\nabla \phi_{K_3},\nabla \phi_{K_1})_K & (\nabla \phi_{K_3},\nabla \phi_{K_2})_K & (\nabla \phi_{K_3},\nabla \phi_{K_3})_K
\end{bmatrix},\end{equation} for all $K\in \mathcal{T}_h$, we end up having all the terms required to perform the summation \eqref{eq:summand}. ($K_j$ is the $j$'th node of triangle $K$) However, we require some sort of knowledge how these so-called "local stiffness matrices" should be summed into the "global stiffness matrix" $S$. For this procedure (often referred to as "assembly of the global stiffness matrix") we define a operator $\mathcal{A}_K: \mathbb{R}^{3\times 3} \rightarrow \mathbb{R}^{N \times N}$ such that
$$(\mathcal{A}_K(X))_{ij} = \begin{cases}
X_{kl}, &\text{if $K_k =i$ and $K_l = j$ for $k,l = 1,2,3$}\\
0,&\text{otherwise.}
\end{cases}$$By using this "assembly operator" we may write
$$S = \sum_{K \in \mathcal{T}_h} \mathcal{A}_K(\hat{S}_K).$$

Let us start coding. Today's language of choice is Mathematica. Unfortunately, the author is not a very experienced Mathematica hacker and therefore the quality of code should be questioned. Any improvements I hear gladly in the commenting section.

We store our mesh node points in a N-by-2 matrix where the first column contains the x-coordinates and the second column the y-coordinates. Thus, let us generate nodes uniformly on the unit square.

mesh[p] =
  Flatten[Outer[{#1,#2}&, Range[0,1,0.1], Range[0,1,0.1]], 1];

Next we need to represent our triangles somehow. For that we use a N-by-3 matrix where each row $j$ contains indices of the $j$'th triangle corresponding to the row indices of the node point matrix. I could have chosen to use the transpose of such matrix instead but I think this is more convenient in the Mathematica-context.

mesh[t] = DelaunayFast[mesh[p]];

The function DelaunayFast generates a Delaunay triangulation of the points in mesh[p]. The function uses the built-in function ListDensityPlot and the implementation I found from StackOverflow reads

Needs["ComputationalGeometry`"];
DelaunayFast[points_] := 
 Module[{pl},
  pl = ListDensityPlot[ArrayPad[points,{{0,0},{0,1}}]];
  Cases[pl, Polygon[a_] :> a, Infinity][[1]]]

In order to visualize the mesh I made an implementation of the function TriMesh which is more or less familiar from the MATLAB environment.

LeaveIfContains[x_, y_] := Select[x, MemberQ[#, y] &];
Tri2Graph[t_] := {#,Union[Flatten[LeaveIfContains[t, #]]]}& /@
   Union[Flatten[t]];
TriMesh[mesh_] := PlanarGraphPlot[mesh[p], Tri2Graph[mesh[t]]];

Now by running

TriMesh[mesh]

we get a picture like

A triangular mesh of the unit square.

For convenience (and to fight Mathematica's annoying [[]]-syntax), we define the following shorthands

mesh[p, k_] := mesh[p][[k]]
mesh[t, k_, l_] := mesh[t][[k, l]]

Then the code to generate all the affine mappings \eqref{eq:affmap} is

A = Table[{mesh[p,mesh[t,k,2]] - mesh[p,mesh[t,k,1]], 
     mesh[p,mesh[t,k,3]] - mesh[p,mesh[t,k,1]]}//Transpose,
     {k, 1, Length[mesh[t]]}];

It is generally a good idea to try to vectorize any interpreted code for performance. Even though this is the first time I'm implementing a solver for Mathematica, I decide to follow the same route I'd use in MATLAB or Numpy: I will write a function which evaluates the elements (i,j) of the local stiffness matrices for all elements at once. To actually evaluate the elements of $\hat{S}_K$ we will transform the integration into the local reference element and use a numerical quadrature. On the reference element the linear basis functions corresponding to the nodes $(0,0), (1,0)$ and $(0,1)$ are $$\begin{aligned} \hat{\phi}_1 &= 1-\hat{x}-\hat{y},\\ \hat{\phi}_2 &= \hat{x}, \\ \hat{\phi}_3 &= \hat{y}. \end{aligned}$$ Exercise: Show that the transformation formula for the gradients is $$\nabla \phi_{K_i} = A_K^{-T} \hat{\nabla} \hat{\phi}_i,$$where $\hat{\nabla}$ denotes the gradient in the local coordinates and $i=1,2,3$. $\blacksquare$

Because the basis functions are linear and defined on a triangle, we may just use a very simple numerical quadrature (which is exact for linear $f$): $$\int_{\hat{K}} f(\hat{x},\hat{y})\,\mathrm{d}x\,\mathrm{d}y \approx \frac 12 f(\frac 13, \frac 13).$$ Thus, the element of the local stiffness matrix is $$(\hat{S}_K)_{ij} = (A_K^{-T} \hat{\nabla} \hat{\phi}_i,A_K^{-T} \hat{\nabla} \hat{\phi}_j)_{\hat{K}} |\det A| = \frac 12\langle A_K^{-T} \hat{\nabla} \hat{\phi}_i,A_K^{-T} \hat{\nabla} \hat{\phi}_j\rangle |\det A|,$$ where $\langle.\rangle$ denotes the standard vector dot product. Now we just may define our local basis as

Phi = {1 - #1 - #2 &, #1 &, #2 &};
DPhi = {{-1., -1.}, {1., 0.}, {0., 1.}};

and then the code to build the elements $(i,j)$ of the local stiffness matrices (for all elements at once) reads

Shat[i_, j_] := 
  0.5 Dot[Inverse[Transpose[#]].DPhi[[i]], 
      Inverse[Transpose[#]].DPhi[[j]]]*Abs[Det[#]] & /@ A;

Similarly for the local load vector, i.e. $$(\hat{f}_K)_i = \frac 12 F(A_K \frac 13 + b_K) \hat{\phi}_i(\frac 13, \frac 13) | \det A_k |,$$ we have

F[x_,y_] := 1.0;
b = Table[mesh[p,mesh[t,k,1]], {k, 1, Length[mesh[t]]}];
fhat[i_] := Phi[[i]][1/3, 1/3]*
    Thread[X[A,b]]/.{X->(0.5 F@@(#1.{1/3,1/3}+#2)*Abs[Det[#1]]&)};

(Sorry for the confusing substitution including X. For some reason Thread didn't work as I expected based on the manual, so I had to come up with a work-around.)

After these function definitions, the global stiffness matrix and load vector assembly is as simple as this:

S = Sum[SparseArray[
    Table[{mesh[t,k,i],mesh[t,k,j]}, {k,1,Length[mesh[t]]}] -> 
    Shat[i,j],{Length[mesh[p]],Length[mesh[p]]}],{i,1,3},{j,1,3}];
f = Sum[SparseArray[
    Table[{mesh[t,k,i],1},{k,1,Length[mesh[t]]}] -> 
    fhat[i],{Length[mesh[p]],1}],{i,1,3}];

Note! This of course assumes that we have set

SetSystemOptions["SparseArrayOptions" ->
    "TreatRepeatedEntries"->Total];

since we get double index pair entries into the SparseArray function call.

We're left with finding the Dirichlet nodes (the boundary nodes) and then solving the system in the inner nodes. To find the boundary nodes I wrote the following short query.

ID = 
  Union[Flatten[Position[mesh[p], #]]] & /@ {{0.,_}, {1.,_},
   {_,0.}, {_, 1.}} // Flatten // Union;

The inner nodes are just a complement of this set:

II = 
  Complement[Range[1, Length[mesh[p]]], ID];

Now we can solve the system by

u = ConstantArray[0., Length[mesh[p]]];
u[[II]] = LinearSolve[S[[II, II]],f[[II]] // Flatten;

and we can visualize the solution by using an implementation of MATLAB's TriSurf I wrote for the occasion:

DrawTriangle[st_, p_, u_] := 
  Graphics3D[Polygon[Append[p[[#]], u[[#]]] & /@ st]];
TriSurf[mesh_, u_] := 
  Show[DrawTriangle[#, mesh[p], u] & /@ mesh[t],
   PlotRange -> All, Axes -> True, BoxRatios -> {1,1,.5}];

By running

TriSurf[mesh,u]

we get the beautiful picture

The solution $u_h$ of the Poisson equation with unit load.

To be continued ...

An inner product defined by a positive definite matrix

Suppose we have positive definite, real and symmetric matrix $D \in \mathbb{R}^{2x2}$. One might be interested in showing the following theorem.

Theorem 1: There exists a constant $C>0$ such that
$$x^T D x \geq C x^T x,$$ holds for every $x \in \mathbb{R}^2$.

This is not as trivial as it seems. To do the trick, one needs to know a few preliminary results. From now on we assume that we're dealing with vectors of $\mathbb{R}^2$ and the matrix $D$ is always $2\times 2$.

Theorem 2: Eigenvalues of a positive definite matrix are positive.

Proof: By definition, positive definiteness means
$$ x^T D x > 0,$$ for every $x$. Then for the eigenvectors $v$ and eigenvalues $\lambda$ of $D$ it holds that
$$ Dv = \lambda v,$$ which can be multiplied from left by $v^T$ to arrive at
$$ v^TDv = \lambda v^T v = \lambda \| v \|^2 > 0,$$ where we used the definition of positive definiteness. Thus, $\lambda > 0$. $\blacksquare$

Theorem 3: A real symmetric matrix $D$ can be decomposed as
$$D = Q \Lambda Q^T,$$ where $Q$ is orthogonal and $\Lambda$ is a diagonal matrix with eigenvalues of $D$ on the diagonal. This is known as the eigendecomposition of the matrix $D$.

Proof: Look up any book on matrix decompositions. $\blacksquare$

Then we may proceed to the original problem.

Proof of Theorem 1: Consider
$$ \frac{x^T D x}{x^T x} = \frac{x^T Q\Lambda Q^T x}{\| x \|^2} = y^T \lambda y,$$ where $y := Q^T \frac{x}{\|x\|}$. The length of $\frac{x}{\|x\|}$ is unity and so is the length of $y$ because $Q^T$ is orthogonal and such a mapping is either a rotation or a mirroring. Thus,
$$y^T \Lambda y = \lambda_1 y_1^2 + \lambda_2 y_2^2 \leq \lambda_{\text{min}} ( y_1^2 + y_2^2) = \lambda_{\text{min}},$$ where $\lambda_{\text{min}} = \min_j \lambda_j, j=1,2$. Multiplying with $x^T x$ gives the result. $\blacksquare$