Deep dive: evaluate state using sparse matrices
Let us write again a weak formulation of the model:
(1)\[\begin{split}\begin{aligned}
\frac{\mathrm{d}}{\mathrm{d}t} \langle y (t), \varphi \rangle_H + \mu_1 a_1 (y(t), \varphi) - \mu_2 \langle f(y(t), q(t)), \varphi \rangle_{H} &= 0 \qquad \forall \varphi \in V, \quad t>0 \\
y(0) &= y_\circ \qquad \text{in } H, \\
\mu_3 a_2 (q(t), \psi) + \mu_4 \langle f(y(t), q(t)), \psi \rangle_{H} - \langle b(t), \psi \rangle_{V'_\circ, V_\circ} &= 0 \qquad \forall \psi \in V_\circ, \quad t>0.
\end{aligned}\end{split}\]
Given approximated spaces \(V^h = \mathrm{span} \, \{ \varphi_1, \dots, \varphi_{N_\mathsf{y}} \} \subset V\) and \(V_\circ^h = \mathrm{span} \, \{ \psi_1, \dots, \psi_{N_\mathsf{y}} \} \subset V_\circ\), the finite element (FE) method is based on the Galerkin projection, i.e. finding approximated solutions \(y^h\) of \(y\) and \(q^h\) of \(q\) of the form
\[\begin{equation*}
y^h (t) = \sum_{i=1}^{N_\mathsf{y}} y^h_i (t) \varphi_i(x) \qquad \text{and} \qquad q^h (t) = \sum_{i=1}^{N_\mathsf{q}} q^h_i (t) \psi_i(x)
\end{equation*}\]
solving the system (1) on the discretized spaces \(V^h\) and \(V_\circ^h\). This is:
(2)\[\begin{split}\begin{aligned}
\frac{\mathrm{d}}{\mathrm{d}t} \langle y^h (t), \varphi_i \rangle_H + \mu_1 a_1 (y^h (t), \varphi_i) - \mu_2 \langle f(y^h(t), q^h(t)), \varphi_i \rangle_{H} &= 0 \qquad i=1,\dots,N_\mathsf{y}, \quad t>0, \\
\langle y^h(0) - y_\circ, \varphi_i \rangle_H &= 0 \qquad i=1,\dots,N_\mathsf{y}, \\
\mu_3 a_2 (q^h (t), \psi_i) + \mu_4 \langle f(y^h(t), q^h(t)), \psi_i \rangle_{H} - \langle b(t), \psi_i \rangle_{V'_\circ, V_\circ} &= 0 \qquad i=1,\dots,N_\mathsf{q}, \quad t>0.
\end{aligned}\end{split}\]
We can define the coordinate vectors
\[\begin{equation*}
\mathbf{y}^h (t) = \left[ \begin{array}{c} y^h_1 (t) \\ \vdots \\ y^h_{N_\mathsf{y}} (t) \end{array} \right] \in \mathbb{R}^{N_\mathsf{y}} \qquad \text{and} \qquad \mathbf{q}^h (t) = \left[ \begin{array}{c} q^h_1 (t) \\ \vdots \\ q^h_{N_\mathsf{q}} (t) \end{array} \right] \in \mathbb{R}^{N_\mathsf{q}}
\end{equation*}\]
At this point, we can approximate (2) in matrix form as:
where:
\[\begin{split}& \left( \mathrm{M}_\mathsf{y} \right)_{i,j} = \langle y^h (t), \varphi_i \rangle_H = \int_\Omega \varphi_j \varphi_i \mathrm{d}x, \\
& \left( \mathrm{A}_\mathsf{y} \right)_{i,j} = a^1_{\boldsymbol\mu} (y^h (t), \varphi_i) = \int_\Omega \kappa_1 \varphi'_j \varphi'_i \mathrm{d}x, \\
& \left( \mathrm{A}_\mathsf{q} \right)_{i,j} = a^2_{\boldsymbol\mu} (y^h (t), \varphi_i) = \int_\Omega \kappa_2 \psi'_j \psi'_i \mathrm{d}x, \\
& \left( \mathbf{y}^h_\circ \right)_i = \langle y_\circ, \varphi_i \rangle_H = \int_\Omega y_\circ \varphi_i \mathrm{d}x, \\
& \left( \mathbf{b}^h (t) \right)_i = \langle b(t), \psi_i \rangle_{V'_\circ, V_\circ},\end{split}\]
and the function \(f\) is evaluated component by component:
\[\left[ f \left( \mathbf{y}^h (t), \mathbf{q}^h (t) \right) \right]_i = f \left( y^h_i (t), q^h_i (t) \right) = \sqrt{y^h_i (t)}\sinh{q^h_i (t)}.\]
To solve (3) we will next discretize in time and then solve the nonlinear problem using Newton’s method. We operate on a uniform time grid \(t_k = k \delta\) for \(k=0, \dots, K\), and \(\delta = T/K\), and we interpret \(\mathbf{y}^{h,k}\) the solution at time \(t_k\) of the fully approximated finite-dimensional system.
In particular, implicit Euler (IE) discretization in time gives:
Let us define the coupled variable
\[\begin{split}\mathbf{z}^{h,k} =
\left(
\begin{array}{c}
\mathbf{y}^{h,k} \\
\mathbf{q}^{h,k}
\end{array}
\right) \in \mathbb{R}^{N_\mathsf{y} + N_\mathsf{q}}\end{split}\]
and the functional
\(\mathrm{F}^k : \mathbb{R}^{N_\mathsf{y} + N_\mathsf{q}} \to \mathbb{R}^{N_\mathsf{y} + N_\mathsf{q}}\), defined as
\[\begin{split}\mathrm{F}^k \left( \mathbf{z}^{h,k} \right) =
\left(
\begin{array}{c}
\mathrm{F}^k_\mathsf{y} \left( \mathbf{z}^{h,k} \right) \\
\mathrm{F}^k_\mathsf{q} \left( \mathbf{z}^{h,k} \right)
\end{array}
\right),\end{split}\]
with
\[\begin{split}& \mathrm{F}^k_\mathsf{y} \left( \mathbf{z}^{h,k} \right) = \left( \mathrm{M}_\mathsf{y} + \mu_1 \delta \mathrm{A}_\mathsf{y} \right) \mathbf{y}^{h,k} - \mu_2 \delta f \left( \mathbf{y}^{h,k}, \mathbf{q}^{h,k} \right) - \mathrm{M}_\mathsf{y} \mathbf{y}^{h,k-1}, \\
& \mathrm{F}^k_\mathsf{q} \left( \mathbf{z}^{h,k} \right) = \mu_3 \mathrm{A}_\mathsf{q} \mathbf{q}^{h,k} + \mu_4 f(\mathbf{y}^{h,k}, \mathbf{q}^{h,k}) - \mathbf{b}^{h,k}.\end{split}\]
Newton’s method for solving \(\mathrm{F}^k (\mathbf{z}^{h,k}) = 0\) at each time step \(k\) works as follow:
Start from \(n=0\) and some initial value \(\mathbf{z}_n\).
Evaluate \(\mathbf{d}\) solution of \(\mathrm{J}_{\mathrm{F}^k} (\mathbf{z}_n) \mathbf{d} = - \mathrm{F}^k(\mathbf{z_n})\).
Generate new point \(\mathbf{z}_{n+1} = \mathbf{z}_n + \mathbf{d}\) and \(n=n+1\).
If \(\mathrm{F}^k(\mathbf{z}_{n})\) is smaller than some tolerance \(\tau\), then we set \(\mathbf{z}^{h,k} := \mathbf{z}_{n}\).
Additionally, at iteration \(k\) we set the previous iterate \(\mathbf{z}^{h,k-1}\) as the initial point \(\mathbf{z}_0\) and we use the variation “damped Newton method” to guarantee monotonicity in the residual decrease.
The overall algorithm to find \(\left \{ \mathbf{y}^{h,k}, \mathbf{q}^{h,k} \right \}_{k=1, \dots, K}\) is:
Evaluate \(\mathbf{y}^{h,0}\) by projecting the intial value \(y_\circ\) onto the solution space.
Evaluate a consistent initial condition for \(\mathbf{q}^{h,0}\) (this step is optional).
Use Newton’s method to evaluate \(\mathbf{z}^{h,k}\) for \(k=1,\dots, K\).