diff --git a/doc/relaxed-micromorphic-continuum.bib b/doc/relaxed-micromorphic-continuum.bib index 30312664f1b9fd5822f07351c4aedeff0171b75e..4123627d14c15ca06a5a1153b63822beea3be967 100644 --- a/doc/relaxed-micromorphic-continuum.bib +++ b/doc/relaxed-micromorphic-continuum.bib @@ -8,3 +8,11 @@ pages = {53--84}, doi = {10.1093/qjmam/hbu027} } + +@Book{hairer_lubich_wanner:2006, + title = {Geometric Numerical Integration---Structure-Preserving Algorithms for Ordinary Differential Equations}, + publisher = {Springer}, + year = {2016}, + author = {Ernst Hairer and Christian Lubich and Gerhard Wanner}, + edition = {Zweite Auflage}, +} diff --git a/doc/relaxed-micromorphic-continuum.tex b/doc/relaxed-micromorphic-continuum.tex index dd1c640b5567447fd71636cd709824eaccf8ec06..6d87ee53f081027e84c511c8eadf97a7d469ae52 100644 --- a/doc/relaxed-micromorphic-continuum.tex +++ b/doc/relaxed-micromorphic-continuum.tex @@ -88,10 +88,9 @@ matrix curl as (\Curl P)_i \colonequals \Big( \frac{\partial P_{i3}}{\partial x_2} - \frac{\partial P_{i2}}{\partial x_3}, \quad - \frac{\partial P_{i3}}{\partial x_2} - \frac{\partial P_{i3}}{\partial x_2}, \quad - \frac{\partial P_{i3}}{\partial x_2} - \frac{\partial P_{i3}}{\partial x_2} \Big). + \frac{\partial P_{i1}}{\partial x_3} - \frac{\partial P_{i3}}{\partial x_1}, \quad + \frac{\partial P_{i2}}{\partial x_1} - \frac{\partial P_{i1}}{\partial x_2} \Big). \end{equation*} -\todosander{Wie würde das Modell in 2d aussehen?} We define the quantities \begin{alignat*}{2} @@ -224,5 +223,138 @@ norm $\tnorm{\cdot}_\mathcal{X}$. \section{The dynamic case} +We now investigate a time-dependent relaxed micromorphic continuum, using the +Hamilton formalism. As shown above, the model is described by a quadratic potential energy, +and we assume the kinetic energy to be quadratic as well. + +\subsection{The implicit midpoint rule for quadratic Hamiltonians} + +Consider a model with quadratic potential energy +\begin{equation*} + U(q) = \frac{1}{2} q^TA q + b^T q +\end{equation*} +and quadratic kinetic energy +\begin{equation*} + T(q, \dot{q}) = T(\dot{q}) = \frac{1}{2} \dot{q}^T M \dot{q}, +\end{equation*} +where $A$ and $M$ are symmetric and positive definite matrices. + +As the kinetic energy is quadratic, the corresponding Hamiltonian +is simply the sum of the two energies~\cite[Chapter~VI.1.2]{hairer_lubich_wanner:2006} +Using that $p = M \dot{q}$ we get +\begin{equation*} + H(p,q) + = + T(\dot{q}(p)) + U(q) + = + \frac{1}{2} p^TM^{-T} p + \frac{1}{2} q^T A q + b^T q. +\end{equation*} + +We want to integrate this system using an implicit midpoint rule. As is +well known, such an integrator is symplectic, second-order accurate~\cite[Theorem~VI.3.5]{hairer_lubich_wanner:2006} +and B-stable~\cite{}. One step of this method is +\begin{equation*} + \begin{pmatrix} p^{k+1} \\ q^{k+1} \end{pmatrix} + = + \begin{pmatrix} p^k \\ q^k \end{pmatrix} + + + \tau J^{-1} \nabla H\Big(\frac{p^k + p^{k+1}}{2}, \frac{q^k + q^{k+1}}{2}\Big), +\end{equation*} +where $\tau$ is the time step size, and $J^{-1}$ is the matrix +\begin{equation*} + J^{-1} + \colonequals + \begin{pmatrix} 0 & -I \\ I & 0 \end{pmatrix}. +\end{equation*} +With +\begin{equation*} + \nabla H(p,q) + = + \begin{pmatrix} + M^{-T} p \\ Aq + b + \end{pmatrix} +\end{equation*} +the integrator becomes +\begin{equation*} + \begin{pmatrix} p^{k+1} \\ q^{k+1} \end{pmatrix} + - \frac{\tau}{2} \begin{pmatrix} 0 & -I \\ I & 0 \end{pmatrix} + \underbrace{ + \begin{pmatrix} + M^{-T} p^{k+1} \\ Aq^{k+1} + \end{pmatrix} + }_{ = \begin{pmatrix} M^{-T} & 0 \\ 0 & A \end{pmatrix} \begin{pmatrix} p^{k+1} \\ q^{k+1} \end{pmatrix} } + = + \begin{pmatrix} p^k \\ q^k \end{pmatrix} + + \frac{\tau}{2} \begin{pmatrix} 0 & -I \\ I & 0 \end{pmatrix} + \begin{pmatrix} + M^{-T} p^k \\ Aq^k + 2b + \end{pmatrix} +\end{equation*} +This can be transformed to +\begin{equation*} + \Bigg[ + \begin{pmatrix} I & 0 \\ 0 & I \end{pmatrix} + - \frac{\tau}{2} + \begin{pmatrix} + 0 & -A \\ M^{-T} & 0 + \end{pmatrix} + \Bigg] + \begin{pmatrix} p^{k+1} \\ q^{k+1} \end{pmatrix} + = + \begin{pmatrix} + p^k - \frac{\tau}{2} A q^k - \tau b\\ + q^k + \frac{\tau}{2} M^{-T} p^k + \end{pmatrix}, +\end{equation*} +and this in turn to +\begin{equation*} + \begin{pmatrix} + I & \frac{\tau}{2} A \\ + - \frac{\tau}{2} M^{-T} & I + \end{pmatrix} + \begin{pmatrix} p^{k+1} \\ q^{k+1} \end{pmatrix} + = + \begin{pmatrix} + p^k - \frac{\tau}{2} A q^k - \tau b\\ + q^k + \frac{\tau}{2} M^{-T} p^k + \end{pmatrix}, +\end{equation*} +Eliminate the lower left block: +\begin{equation*} + \begin{pmatrix} + I & \frac{\tau}{2} A \\ + 0 & I + \frac{\tau^2}{4} M^{-T}A + \end{pmatrix} + \begin{pmatrix} p^{k+1} \\ q^{k+1} \end{pmatrix} + = + \begin{pmatrix} + p^k - \frac{\tau}{2} A q^k - \tau b\\ + q^k + \frac{\tau}{2} M^{-T} p^k + \frac{\tau}{2} M^{-T}(p^k - \frac{\tau}{2}Aq^k - \tau b) + \end{pmatrix}. +\end{equation*} + +This can be solved by backward substitution. In algorithmic form this reads: +\begin{enumerate} + \item Solve + \begin{equation*} + \Big(M^T + \frac{\tau^2}{4} A \Big) q^{k+1} + = + M^Tq^k + \tau p^k - \frac{\tau^2}{4}Aq^k - \frac{\tau^2}{2}b + \end{equation*} + + \item Set + \begin{equation*} + p^{k+1} + = + p^k - \frac{\tau}{2} Aq^k - 2b - \frac{\tau}{2}Aq^{k+1}. + \end{equation*} + +\end{enumerate} + + + + + + \printbibliography \end{document}