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.
Tilaa:
Blogitekstit (Atom)