{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e48786b3-d807-411a-b390-4e12d30d8ec0",
   "metadata": {},
   "source": [
    "# Deep dive: evaluate state using sparse matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6651c125-e5a1-4015-b0e9-d9c0460874dd",
   "metadata": {},
   "source": [
    "Let us write again a weak formulation of the model:\n",
    "```{math}\n",
    ":label: weak-model\n",
    "\\begin{aligned}\n",
    "    \\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 \\\\\n",
    "    y(0) &= y_\\circ \\qquad \\text{in } H, \\\\\n",
    "    \\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.\n",
    "\\end{aligned}\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e1874f0-e189-4efc-9c41-b2730e772ba4",
   "metadata": {},
   "source": [
    "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\n",
    "\\begin{equation*}\n",
    "    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)\n",
    "\\end{equation*}\n",
    "solving the system [](#weak-model) on the discretized spaces $V^h$ and $V_\\circ^h$. This is:\n",
    "```{math}\n",
    ":label: galerkin-projection\n",
    "\\begin{aligned}\n",
    "    \\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, \\\\\n",
    "    \\langle y^h(0) - y_\\circ, \\varphi_i \\rangle_H &= 0 \\qquad i=1,\\dots,N_\\mathsf{y}, \\\\\n",
    "    \\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.\n",
    "\\end{aligned}\n",
    "```\n",
    "We can define the coordinate vectors\n",
    "\\begin{equation*}\n",
    "    \\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}}\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "845bfb99-6d60-4123-8acd-b98f781c187a",
   "metadata": {},
   "source": [
    "At this point, we can approximate [](#galerkin-projection) in matrix form as:\n",
    "```{math}\n",
    ":label: galerkin-projection-matrix-form\n",
    "\\begin{aligned}\n",
    "    \\mathrm{M}_\\mathsf{y} \\dot{\\mathbf{y}}^h (t) + \\mu_1 \\mathrm{A}_\\mathsf{y} \\mathbf{y}^h (t) - \\mu_2 f \\left( \\mathbf{y}^h (t), \\mathbf{q}^h (t) \\right) &= 0 \\qquad \\text{for } t>0, \\\\\n",
    "    \\mathrm{M}_\\mathsf{y} \\mathbf{y}^h (0) &= \\mathbf{y}^h_\\circ, \\\\\n",
    "    \\mu_3 \\mathrm{A}_\\mathsf{q} \\mathbf{q}^h (t) + \\mu_4 f(\\mathbf{y}^h (t), \\mathbf{q}^h (t)) &= \\mathbf{b}^h (t) \\qquad \\text{for } \\quad t>0,\n",
    "\\end{aligned}\n",
    "```\n",
    "where:\n",
    "```{math}\n",
    "& \\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, \\\\\n",
    "& \\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, \\\\\n",
    "& \\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, \\\\\n",
    "& \\left( \\mathbf{y}^h_\\circ \\right)_i = \\langle y_\\circ, \\varphi_i \\rangle_H = \\int_\\Omega y_\\circ \\varphi_i \\mathrm{d}x, \\\\\n",
    "& \\left( \\mathbf{b}^h (t) \\right)_i = \\langle b(t), \\psi_i \\rangle_{V'_\\circ, V_\\circ},\n",
    "```\n",
    "and the function $f$ is evaluated component by component:\n",
    "```{math}\n",
    "\\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)}.\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5ae7bc6-de89-4acf-b56e-d8ca341b9f0f",
   "metadata": {},
   "source": [
    "To solve [](#galerkin-projection-matrix-form) 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."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5518600-7fed-49d7-8387-6f4eb4d97366",
   "metadata": {},
   "source": [
    "In particular, implicit Euler (IE) discretization in time gives:\n",
    "```{math}\n",
    ":label: galerkin-projection-matrix-form-ie\n",
    "\\begin{aligned}\n",
    "    \\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} \\qquad \\text{for } k = 1, \\dots, K, \\\\\n",
    "    \\mathrm{M}_\\mathsf{y} \\mathbf{y}^{h,0} &= \\mathbf{y}^h_\\circ, \\\\\n",
    "    \\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} \\qquad \\text{for } k = 1, \\dots, K,\n",
    "\\end{aligned}\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce3f0e7a-ccaf-48ce-a4ad-4a59891f34cf",
   "metadata": {},
   "source": [
    "Let us define the coupled variable\n",
    "```{math}\n",
    "\\mathbf{z}^{h,k} =\n",
    "\\left(\n",
    "  \\begin{array}{c}\n",
    "    \\mathbf{y}^{h,k} \\\\\n",
    "    \\mathbf{q}^{h,k}\n",
    "  \\end{array}\n",
    "\\right) \\in \\mathbb{R}^{N_\\mathsf{y} + N_\\mathsf{q}}\n",
    "```\n",
    "and the functional\n",
    "$\\mathrm{F}^k : \\mathbb{R}^{N_\\mathsf{y} + N_\\mathsf{q}} \\to \\mathbb{R}^{N_\\mathsf{y} + N_\\mathsf{q}}$, defined as\n",
    "```{math}\n",
    "\\mathrm{F}^k \\left( \\mathbf{z}^{h,k} \\right) =\n",
    "\\left(\n",
    "  \\begin{array}{c}\n",
    "    \\mathrm{F}^k_\\mathsf{y} \\left( \\mathbf{z}^{h,k} \\right) \\\\\n",
    "    \\mathrm{F}^k_\\mathsf{q} \\left( \\mathbf{z}^{h,k} \\right)\n",
    "  \\end{array}\n",
    "\\right),\n",
    "```\n",
    "with\n",
    "```{math}\n",
    "& \\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}, \\\\\n",
    "& \\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}.\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52310d64-0fdc-41b9-aa48-3f9e73b2a7c4",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "Newton's method for solving $\\mathrm{F}^k (\\mathbf{z}^{h,k}) = 0$ at each time step $k$ works as follow:\n",
    "- Start from $n=0$ and some initial value $\\mathbf{z}_n$.\n",
    "- Evaluate $\\mathbf{d}$ solution of $\\mathrm{J}_{\\mathrm{F}^k} (\\mathbf{z}_n) \\mathbf{d} = - \\mathrm{F}^k(\\mathbf{z_n})$.\n",
    "- Generate new point $\\mathbf{z}_{n+1} = \\mathbf{z}_n + \\mathbf{d}$ and $n=n+1$.\n",
    "- If $\\mathrm{F}^k(\\mathbf{z}_{n})$ is smaller than some tolerance $\\tau$, then we set $\\mathbf{z}^{h,k} := \\mathbf{z}_{n}$.\n",
    "\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."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9de86278-1d31-409b-889a-2bdd779e040f",
   "metadata": {},
   "source": [
    "The overall algorithm to find $\\left \\{ \\mathbf{y}^{h,k}, \\mathbf{q}^{h,k} \\right \\}_{k=1, \\dots, K}$ is:\n",
    "- Evaluate $\\mathbf{y}^{h,0}$ by projecting the intial value $y_\\circ$ onto the solution space.\n",
    "- Evaluate a consistent initial condition for $\\mathbf{q}^{h,0}$ (this step is optional).\n",
    "- Use Newton's method to evaluate $\\mathbf{z}^{h,k}$ for $k=1,\\dots, K$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dcf15d3b-cb85-4445-bd5e-239a6d00a5c3",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a5f1e1d-baee-44c6-8b7a-4a78101250fe",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86f65ce9-1a86-4ab2-b10a-243dbbeea46f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c291adfe-3505-4292-b22d-56ab6ecd30c2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e49e73d-b61b-42c1-b08a-4a926be6a6a4",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}