---
title: Eigenanalysis
---

# Eigen Analysis for the Cantilever Bar Example

We aim to analyze and explain the dynamic behavior through eigen analysis, focusing on the cantilever bar example, as we have a clear understanding of how the leading eigenvector should look like. To figure out why the baseline progressive version damps even more motion compared to the direct coarse simulation, we compute and compare the eigenvalues corresponding to the same eigenvector, where larger eigenvalues suggeest increased numerical stiffening.

### :stars: Theory
Suppose we have a particle system with potential $V(\mathbf{q})$, so that the equations of motion are given by
$$M\ddot{\mathbf{q}} = -[dV(\mathbf{q})]^T.$$

Let us now assume that $\mathbf{q_0}$ is a stable equilibrium of the system. Suppose we want to simulate the system, starting from a configuration near $\mathbf{q_0}$. Instead of using $\mathbf{q}$ as our degrees of freedom, let us introduce new degrees of freedom $\mathbf{d}$, _displacements away from_ $\mathbf{q_0}$, related to the original coordinates by
$$\mathbf{q} = \mathbf{q_0} + \mathbf{d} \quad \Rightarrow \ \ M\ddot{\mathbf{d}} = -[dV(\mathbf{q_0 + d})]^T.$$

By Taylor expansion and after dropping off the higher-order terms, we get
\begin{equation}
M\ddot{\mathbf{d}} \approx -HV(\mathbf{q_0})\mathbf{d}.
\end{equation}

This is an ODE with analytic solutions which boils down to solving the generalized eigenvalue problem
$$HV(\mathbf{q_0})\mathbf{v_i} = \lambda_i M\mathbf{v_i},$$

where $\lambda_i$ and $\mathbf{v_i}$ be the generalized eigenvalues and eigenvectors of the hessian $HV(\mathbf{q_0})$ with respect to the mass matrix $M$.

For our cantilever bar example, we compute the hessian $H$ at the rest shape (must be **SPD**) where $H$ is the sum of the membrane hessian and bending hessian.

<p>
<img style="float: right;" src="https://hackmd.io/_uploads/r1eMfzu3h.png" width="45%">
</p>

Red dots on the left show fixed dirichlet boundary constraints. Because of these, when computing the modes, we also need to enforce an additional equality constraint $S^T \mathbf{v_i} = 0$, where $S$ is the selection matrix. This can be achieved through a constraint-augmented generalized eigenvalue problem:

\begin{align}
\begin{bmatrix}
\boldsymbol H & \boldsymbol S \\
\boldsymbol S^T & \boldsymbol 0
\end{bmatrix}
\begin{bmatrix}
\boldsymbol{v} \\
\boldsymbol{\mu}
\end{bmatrix}
=
\boldsymbol{\Lambda}
\begin{bmatrix}
\boldsymbol v^T & \boldsymbol \mu ^T
\end{bmatrix}
\begin{bmatrix}
\boldsymbol M & \boldsymbol 0 \\
\boldsymbol 0 & \boldsymbol 0
\end{bmatrix}.
\end{align}

### :stars: Results
In this experiment, we compare the following four different versions of the same setup including material parameters, initial conditions etc. and compare their eigenvalues of the leading eigenmode.

+ V0. Direct fine level simulation
+ Both the hessian and mass matrix are sampled from the fine level
+ $H = H_f$ and $M = M_f$
+ V1. Direct coarse level simulation
+ Both the hessian and mass matrix are sampled from the coarse level
+ $H = H_c$ and $M = M_c$
+ V2. Baseline progressive simulation with the coarse-level mass matrix
+ $P$ is the prolongation matrix
+ $H = P^T * H_f * P$ and $M = M_c$
+ V3. Baseline progressive simulation with the sandwiched fine-level mass matrix
+ $P$ is the prolongation matrix
+ $H = P^T * H_f * P$ and $M = P^T * M_f * P$

As an example, here are the first eight leading eigenmodes for the baseline progressive simulation weighted by the coarse-level mass matrix.
<p>
<img style="float: center;" src="https://hackmd.io/_uploads/r1BXolO3h.png" width="49%">
<img style="float: center;" src="https://hackmd.io/_uploads/HJlHQjeu23.png" width="49%">
<img style="float: center;" src="https://hackmd.io/_uploads/HJHQoedn3.png" width="49%">
<img style="float: center;" src="https://hackmd.io/_uploads/SyxHXogO2h.png" width="49%">
<img style="float: center;" src="https://hackmd.io/_uploads/B1Bmogd23.png" width="49%">
<img style="float: center;" src="https://hackmd.io/_uploads/HJBXsxO22.png" width="49%">
<img style="float: center;" src="https://hackmd.io/_uploads/B1H7jldn3.png" width="49%">
<img style="float: center;" src="https://hackmd.io/_uploads/BJHXsx_nn.png" width="49%">
</p>

And the first eigenmode basically remains consistent across four different versions of simulations:

<p>
<img style="float: center;" src="https://hackmd.io/_uploads/H1gYVM_22.png" width="24%">
<img style="float: center;" src="https://hackmd.io/_uploads/SyOxQzdhh.png" width="24%">
<img style="float: center;" src="https://hackmd.io/_uploads/r1BXolO3h.png" width="24%">
<img style="float: center;" src="https://hackmd.io/_uploads/H1FmBz_3h.png" width="24%">
</p>

The corresponding eigenvalues are shown as below:
|Version |Direct Fine|Direct Coarse|PD Coarse|PD Sandwiched Fine|
|-----|--------|--------|--------|--------|
|$\lambda_0$|0.003405|0.0043165|0.013854|0.013859|

#### Unconstrained version:
|Version |Direct Fine|Direct Coarse|PD Coarse|PD Sandwiched Fine|
|-----|--------|--------|--------|--------|
|$\lambda_0$|0.12897|0.16973|0.56837|0.57052|

which confirms the reason why the motion of the baseline progressive dynamic simulation seems even more damped out (stiffer).

### :stars: Opportunities
+ Rigorous eigen analysis of the subspace parameterized equation of motion:
$$M\ddot{P}(q_c) = -[dV(P(q_c))]^T.$$
where the definition of $M$ is up to our choice, e.g. increase the weight.
+ Alternatively, if we look at its variational form ,i.e., the incremental potential, the definition of the inertia term is also something we could further play with.
$$
x^{t+1} = \underset{x}{\operatorname{argmin}} \frac{1}{2 h^2} \|x - x^t \|^2_{M} + E(x, \bar x, u^{t+1}).
$$
+ Look into why the PD Sandwiched Fine version doesn't help.