keskiviikko 10. joulukuuta 2014

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?

lauantai 6. joulukuuta 2014

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.