An expression for the derivative of the determinant function is useful occasionally in mechanics (e.g. when solving hyperelastic models by the FEM). Derivation of the identity is, however, pretty annoying. Fortunately there is a trick you can employ:
$$\det(A+\epsilon B)=\det(A)\det(I+\epsilon A^{-1}B).$$Let $C=A^{-1}B$. Then
$$\det(I+\epsilon C)=\exp\log\det(I+\epsilon C) = \exp \mathrm{tr} \log(I+\epsilon C).$$Using the series definition of matrix logarithm
$$\log(I+\epsilon C) = \sum_{m=1}^\infty \frac{(-1)^{m+1}(\epsilon C)^m}{m}.$$Thus,
$$\frac{\mathrm{d}}{\mathrm{d} \epsilon} \det(I+\epsilon C) = \mathrm{tr}(\sum_{m=1}^\infty (-1)^{m+1}\epsilon^{m-1} C^m) \exp \mathrm{tr} \log(I+\epsilon C)$$implying
$$D_G(\det(A))B = \det(A) \mathrm{tr}(A^{-1}B) \det(I) =\det(A) \mathrm{tr}(A^{-1}B).$$What we have obtained is the Gateaux derivative of $\det(A)$ at $A$ into direction $B$. In this case the operator $D_G(\det(A))$ is linear (to see this, write $\mathrm{tr}(A^{-1}B)=A^{-T}:B$) and, hence, we have also obtained the Frechet derivative:
$$D_F(\det(A)) = \det(A)A^{-T}$$
Note that $\log(\det(A)) = \mathrm{tr}(\log(A))$ follows by considering the eigenvalues of $A$ (determinant is the product of eigenvalues) and $\log(A)$ (trace is the sum of eigenvalues).
lauantai 7. marraskuuta 2015
keskiviikko 10. kesäkuuta 2015
Notes to myself on stress tensors
Cauchy stress is the true stress on a deformed configuration. It operates on normal vectors on deformed solid and gives traction on deformed solid.
First Piola-Kirchoff tensor operates to normal vectors on undeformed solid and gives traction on deformed solid. (This is also 3D generalization of "engineering stress" and unsymmetric in general.)
Second Piola-Kirchoff tensor operates to normal vectors on undeformed solid and gives the "pullback" of traction on undeformed solid. ("pullback" of traction = traction inverse mapped through deformation.)
First Piola-Kirchoff tensor operates to normal vectors on undeformed solid and gives traction on deformed solid. (This is also 3D generalization of "engineering stress" and unsymmetric in general.)
Second Piola-Kirchoff tensor operates to normal vectors on undeformed solid and gives the "pullback" of traction on undeformed solid. ("pullback" of traction = traction inverse mapped through deformation.)
maanantai 1. kesäkuuta 2015
VIM commands that I found useful
ctrl-d
Autocompletion in normal mode commands. Surprisingly useful with e.g. :help q<ctrl-d>.
q<register><commands>q
Record a macro to <register> (can be any alphabet in a-z). <commands> can include multiple commands.
<times>@<register>
Execute macro <times> times from <register>.
I in block visual mode
Insert text to multiple lines.
c in block visual mode
Change text in multiple lines.
ci{ or ci[ or ci(
"Change inside <brackettype>". Also works with vi{, i.e. "Visual selection inside <brackettype>".
zf (in visual mode) and za (in insert mode)
Create a code fold and open/close the fold. Also :mkview can be used to save the folds and :loadview to load previously saved folds.
ci{ or ci[ or ci(
"Change inside <brackettype>". Also works with vi{, i.e. "Visual selection inside <brackettype>".
zf (in visual mode) and za (in insert mode)
Create a code fold and open/close the fold. Also :mkview can be used to save the folds and :loadview to load previously saved folds.
lauantai 7. maaliskuuta 2015
Newton's method and nonlinear finite elements
Newton-Galerkin method, the most carefully guarded secret of finite element business. All examples I have seen so far feel like 'special tricks' and therefore I want to cover it once-and-for-all thoroughly, even just to remind myself.
Let us first consider searching an approximation to the zeros of a function $f : \mathbb{R} \rightarrow \mathbb{R}$ (i.e. standard Newton's method) because this is very analogous to what we do when solving zeros of functionals. Thus, let us solve
$$f(x)=0.$$ Suppose that we have a sufficiently good approximation $x_0$ to $x$. Then we attempt to compute the correct perturbation $\Delta x$ satisfying
$$x = x_0 + \Delta x.$$ Then using Taylor's expansion we get
$$0 = f(x) = f(x_0 + \Delta x) = f(x_0)+f'(x_0) \Delta x + \mathcal{O}(\Delta x^2).$$ Neglecting the higher order terms and solving for $\Delta x$ gives the approximate equation
$$\Delta x \approx - \frac{f(x_0)}{f'(x_0)}.$$ Then we have a new approximation to $x$:
$$x \approx x_0 - \frac{f(x_0)}{f'(x_0)},$$ which you might remember from high school (back then we used some stupid geometrical argument).
Similarly, let us consider a non-linear PDE for which the weak formulation is: find $u \in V$ such that
$$a(u;v)=l(v),\quad \forall v \in V,$$ where $a(u;v)$ is linear in $v$, non-linear in $u$, and $l(v)$ is linear. Here $V$ is some function space with the proper boundary conditions
Example. A weak form like this could, for example, arise from the strong formulation
$$-\nabla \cdot (u^3 \nabla u) = f,\quad \text{in $\Omega$,}$$ with the boundary condition $u|_{\partial \Omega}=0$. In this case
$$a(u;v)=\int_\Omega u^3 \nabla u \cdot \nabla v \,\mathrm{d}x$$ and the usual
$$l(v)=\int_\Omega f v\,\mathrm{d}x.$$
Anyways, we again assume that we have an approximation $u_0$ to $u$ and search for the correction $\delta u$ satisfying
$$u = u_0 + \delta u.$$ Then a similar Taylor's expansion as before gives
$$0 = a(u_0+\delta u;v) - l(v) = a(u_0;v)+D_u a(u_0;v) \delta u + \mathcal{O}(\delta u^2) - l(v),$$ where $D_u a(u_0;v) \delta u$ is the derivative of $a$ with respect to the first argument into the direction $\delta u$ evaluated at $(u_0,v)$. Rearranging and neglecting the HOT's gives
$$D_u a(u_0;v) \delta u \approx l(v)-a(u_0;v).$$ Note that the derivative $D_u a(u_0;v)$ is a linear operator acting on $\delta u$ and therefore in general we must discretize the previous equation before solving for $\delta u$. However, first we should compute an explicit representation for $D_u a(u_0;v) \delta u$ and this depends on how $a$ is defined. An example follows.
Example. Suppose
$$a(u;v)=\int_\Omega u^3 \nabla u \cdot \nabla v\,\mathrm{d}x$$ as in the previous example. Then by the definition of the directional derivative
\begin{align*}
D_u a(u_0;v) \delta u &= \frac{\partial a(u_0+t \delta u;v)}{\partial t}\Big|_{t=0}\\
&= \frac{\partial}{\partial t} \int_\Omega (u_0 + t \delta u)^3 \nabla(u+t \delta u) \cdot \nabla v\,\mathrm{d}x \Big|_{t=0} \\
&= \int_\Omega 3 u_0^2 \delta u \nabla u_0 \cdot \nabla v\,\mathrm{d}x + \int_\Omega u_0^3 \nabla \delta u \cdot \nabla v\,\mathrm{d}x.
\end{align*}
Now one step of the Newton-Galerkin iteration reads: find $\delta u \in V$ such that
$$D_u a(u_0;v) \delta u = l(v)-a(u_0;v)$$ holds for every $v \in V$.
Let us first consider searching an approximation to the zeros of a function $f : \mathbb{R} \rightarrow \mathbb{R}$ (i.e. standard Newton's method) because this is very analogous to what we do when solving zeros of functionals. Thus, let us solve
$$f(x)=0.$$ Suppose that we have a sufficiently good approximation $x_0$ to $x$. Then we attempt to compute the correct perturbation $\Delta x$ satisfying
$$x = x_0 + \Delta x.$$ Then using Taylor's expansion we get
$$0 = f(x) = f(x_0 + \Delta x) = f(x_0)+f'(x_0) \Delta x + \mathcal{O}(\Delta x^2).$$ Neglecting the higher order terms and solving for $\Delta x$ gives the approximate equation
$$\Delta x \approx - \frac{f(x_0)}{f'(x_0)}.$$ Then we have a new approximation to $x$:
$$x \approx x_0 - \frac{f(x_0)}{f'(x_0)},$$ which you might remember from high school (back then we used some stupid geometrical argument).
Similarly, let us consider a non-linear PDE for which the weak formulation is: find $u \in V$ such that
$$a(u;v)=l(v),\quad \forall v \in V,$$ where $a(u;v)$ is linear in $v$, non-linear in $u$, and $l(v)$ is linear. Here $V$ is some function space with the proper boundary conditions
Example. A weak form like this could, for example, arise from the strong formulation
$$-\nabla \cdot (u^3 \nabla u) = f,\quad \text{in $\Omega$,}$$ with the boundary condition $u|_{\partial \Omega}=0$. In this case
$$a(u;v)=\int_\Omega u^3 \nabla u \cdot \nabla v \,\mathrm{d}x$$ and the usual
$$l(v)=\int_\Omega f v\,\mathrm{d}x.$$
Anyways, we again assume that we have an approximation $u_0$ to $u$ and search for the correction $\delta u$ satisfying
$$u = u_0 + \delta u.$$ Then a similar Taylor's expansion as before gives
$$0 = a(u_0+\delta u;v) - l(v) = a(u_0;v)+D_u a(u_0;v) \delta u + \mathcal{O}(\delta u^2) - l(v),$$ where $D_u a(u_0;v) \delta u$ is the derivative of $a$ with respect to the first argument into the direction $\delta u$ evaluated at $(u_0,v)$. Rearranging and neglecting the HOT's gives
$$D_u a(u_0;v) \delta u \approx l(v)-a(u_0;v).$$ Note that the derivative $D_u a(u_0;v)$ is a linear operator acting on $\delta u$ and therefore in general we must discretize the previous equation before solving for $\delta u$. However, first we should compute an explicit representation for $D_u a(u_0;v) \delta u$ and this depends on how $a$ is defined. An example follows.
Example. Suppose
$$a(u;v)=\int_\Omega u^3 \nabla u \cdot \nabla v\,\mathrm{d}x$$ as in the previous example. Then by the definition of the directional derivative
\begin{align*}
D_u a(u_0;v) \delta u &= \frac{\partial a(u_0+t \delta u;v)}{\partial t}\Big|_{t=0}\\
&= \frac{\partial}{\partial t} \int_\Omega (u_0 + t \delta u)^3 \nabla(u+t \delta u) \cdot \nabla v\,\mathrm{d}x \Big|_{t=0} \\
&= \int_\Omega 3 u_0^2 \delta u \nabla u_0 \cdot \nabla v\,\mathrm{d}x + \int_\Omega u_0^3 \nabla \delta u \cdot \nabla v\,\mathrm{d}x.
\end{align*}
Now one step of the Newton-Galerkin iteration reads: find $\delta u \in V$ such that
$$D_u a(u_0;v) \delta u = l(v)-a(u_0;v)$$ holds for every $v \in V$.
maanantai 23. helmikuuta 2015
Matlab EPS-export fails on dotted lines
I continue my set of tips on drawing bearable figures using Matlab. Have you ever tried to increase the line width of a double dotted (with linespec ':') or dot-dash (with linespec '-.') line? Does it look like crap after exporting to *.eps?
There is a simple solution to the problem: open up the *.eps-file and search for
By increasing the first decimal number 0.5 up to e.g. 1 you should get much better results.
There is a simple solution to the problem: open up the *.eps-file and search for
- "/DO" for ':'
- "/DD" for '-.'
- /DO { [0.5 dpi2point mul 4 dpi2point mul] 0 setdash } bdef
- /DD { [0.5 dpi2point mul 4 dpi2point mul 6 dpi2point mul 4
By increasing the first decimal number 0.5 up to e.g. 1 you should get much better results.
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:
You may compare them to the solutions given at Wikipedia. Can you find them all?
$$-\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:
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.
$$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.
Tilaa:
Blogitekstit (Atom)