7- femsys.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Variational forms for systems of PDEs\n",
    "<div id=\"ch:femsys\"></div>\n",
    "\n",
    "Many mathematical models involve $m+1$ unknown functions\n",
    "governed by a system of $m+1$ differential equations. In abstract form\n",
    "we may denote the unknowns by $u^{(0)},\\ldots,\n",
    "u^{(m)}$ and write the governing equations as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "\\mathcal{L}_0(u^{(0)},\\ldots,u^{(m)}) &= 0,\\\\\n",
    "&\\vdots\\\\\n",
    "\\mathcal{L}_{m}(u^{(0)},\\ldots,u^{(m)}) &= 0,\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $\\mathcal{L}_i$ is some differential operator defining differential\n",
    "equation number $i$.\n",
    "\n",
    "# Variational forms\n",
    "<div id=\"fem:sys:vform\"></div>\n",
    "\n",
    "There are basically two ways of formulating a variational form\n",
    "for a system of differential equations. The first method treats\n",
    "each equation independently as a scalar equation, while the other\n",
    "method views the total system as a vector equation with a vector function\n",
    "as unknown.\n",
    "\n",
    "## Sequence of scalar PDEs formulation\n",
    "\n",
    "Let us start with the approach that treats one equation at a time.\n",
    "We multiply equation number $i$ by\n",
    "some test function $v^{(i)}\\in V^{(i)}$ and integrate over the domain:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:sys:vform:1by1a\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\int_\\Omega \\mathcal{L}^{(0)}(u^{(0)},\\ldots,u^{(m)}) v^{(0)}{\\, \\mathrm{d}x} = 0,\n",
    "\\label{fem:sys:vform:1by1a} \\tag{1}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto1\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "\\vdots\n",
    "\\label{_auto1} \\tag{2}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:sys:vform:1by1b\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "\\int_\\Omega \\mathcal{L}^{(m)}(u^{(0)},\\ldots,u^{(m)}) v^{(m)}{\\, \\mathrm{d}x} = 0\n",
    "\\label{fem:sys:vform:1by1b} \\tag{3}\n",
    "{\\thinspace .}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Terms with second-order derivatives may be integrated by parts, with\n",
    "Neumann conditions inserted in boundary integrals.\n",
    "Let"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "V^{(i)} = \\hbox{span}\\{{\\psi}_0^{(i)},\\ldots,{\\psi}_{N_i}^{(i)}\\},\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "such that"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u^{(i)} = B^{(i)}(\\boldsymbol{x}) + \\sum_{j=0}^{N_i} c_j^{(i)} {\\psi}_j^{(i)}(\\boldsymbol{x}),\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $B^{(i)}$ is a boundary function to handle nonzero Dirichlet conditions.\n",
    "Observe that different unknowns may live in different spaces with different\n",
    "basis functions and numbers of degrees of freedom.\n",
    "\n",
    "From the $m$ equations in the variational forms we can derive\n",
    "$m$ coupled systems of algebraic equations for the\n",
    "$\\Pi_{i=0}^{m} N_i$ unknown coefficients $c_j^{(i)}$, $j=0,\\ldots,N_i$,\n",
    "$i=0,\\ldots,m$.\n",
    "\n",
    "## Vector PDE formulation\n",
    "\n",
    "The alternative method for deriving a variational form for a system of\n",
    "differential equations introduces a vector of unknown functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\boldsymbol{u} = (u^{(0)},\\ldots,u^{(m)}),\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "a vector of test functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\boldsymbol{v} = (v^{(0)},\\ldots,v^{(m)}),\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "with"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\boldsymbol{u}, \\boldsymbol{v} \\in  \\boldsymbol{V} = V^{(0)}\\times \\cdots \\times V^{(m)}\n",
    "{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With nonzero Dirichlet conditions, we have a vector\n",
    "$\\boldsymbol{B} = (B^{(0)},\\ldots,B^{(m)})$ with boundary functions and then\n",
    "it is $\\boldsymbol{u} - \\boldsymbol{B}$ that lies in $\\boldsymbol{V}$, not $\\boldsymbol{u}$ itself.\n",
    "\n",
    "The governing system of differential equations is written"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\boldsymbol{\\mathcal{L}}(\\boldsymbol{u} ) = 0,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\boldsymbol{\\mathcal{L}}(\\boldsymbol{u} ) = (\\mathcal{L}^{(0)}(\\boldsymbol{u}),\\ldots, \\mathcal{L}^{(m)}(\\boldsymbol{u}))\n",
    "{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variational form is derived by taking the inner product of\n",
    "the vector of equations and the test function vector:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:sys:vform:inner\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\int_\\Omega \\boldsymbol{\\mathcal{L}}(\\boldsymbol{u} )\\cdot\\boldsymbol{v} = 0\\quad\\forall\\boldsymbol{v}\\in\\boldsymbol{V}{\\thinspace .}\n",
    "\\label{fem:sys:vform:inner} \\tag{4}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observe that ([4](#fem:sys:vform:inner)) is one scalar equation. To derive\n",
    "systems of algebraic equations for the unknown coefficients in the\n",
    "expansions of the unknown functions, one chooses $m$ linearly\n",
    "independent $\\boldsymbol{v}$ vectors to generate $m$ independent variational forms\n",
    "from ([4](#fem:sys:vform:inner)).  The particular choice $\\boldsymbol{v} =\n",
    "(v^{(0)},0,\\ldots,0)$ recovers ([1](#fem:sys:vform:1by1a)), $\\boldsymbol{v} =\n",
    "(0,\\ldots,0,v^{(m)})$ recovers ([3](#fem:sys:vform:1by1b)), and $\\boldsymbol{v} =\n",
    "(0,\\ldots,0,v^{(i)},0,\\ldots,0)$ recovers the variational form number\n",
    "$i$, $\\int_\\Omega \\mathcal{L}^{(i)} v^{(i)}{\\, \\mathrm{d}x} =0$, in\n",
    "([1](#fem:sys:vform:1by1a))-([3](#fem:sys:vform:1by1b)).\n",
    "\n",
    "# A worked example\n",
    "<div id=\"fem:sys:uT:ex\"></div>\n",
    "\n",
    "We now consider a specific system of two partial differential equations\n",
    "in two space dimensions:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:sys:wT:ex:weq\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\mu \\nabla^2 w = -\\beta,\n",
    "\\label{fem:sys:wT:ex:weq} \\tag{5}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:sys:wT:ex:Teq\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "\\kappa\\nabla^2 T = - \\mu ||\\nabla w||^2\n",
    "{\\thinspace .}\n",
    "\\label{fem:sys:wT:ex:Teq} \\tag{6}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The unknown functions $w(x,y)$ and $T(x,y)$ are defined in a domain $\\Omega$,\n",
    "while $\\mu$, $\\beta$,\n",
    "and $\\kappa$ are given constants. The norm in\n",
    "([6](#fem:sys:wT:ex:Teq)) is the standard Euclidean norm:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "||\\nabla w||^2 = \\nabla w\\cdot\\nabla w = w_x^2 + w_y^2\n",
    "{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The boundary conditions associated with\n",
    "([5](#fem:sys:wT:ex:weq))-([6](#fem:sys:wT:ex:Teq)) are $w=0$ on\n",
    "$\\partial\\Omega$ and $T=T_0$ on $\\partial\\Omega$.\n",
    "Each of the equations ([5](#fem:sys:wT:ex:weq)) and ([6](#fem:sys:wT:ex:Teq))\n",
    "needs one condition at each point on the boundary.\n",
    "\n",
    "The system ([5](#fem:sys:wT:ex:weq))-([6](#fem:sys:wT:ex:Teq)) arises\n",
    "from fluid flow in a straight pipe, with the $z$ axis in the direction\n",
    "of the pipe. The domain $\\Omega$ is a cross section of the pipe, $w$\n",
    "is the velocity in the $z$ direction, $\\mu$\n",
    "is the viscosity of the fluid, $\\beta$ is the pressure gradient along\n",
    "the pipe, $T$ is the temperature,\n",
    "and $\\kappa$ is the heat conduction coefficient of the\n",
    "fluid. The equation ([5](#fem:sys:wT:ex:weq)) comes from the Navier-Stokes\n",
    "equations, and ([6](#fem:sys:wT:ex:Teq)) follows from the energy equation.\n",
    "The term $- \\mu ||\\nabla w||^2$ models heating of the fluid\n",
    "due to internal friction.\n",
    "\n",
    "Observe that the system ([5](#fem:sys:wT:ex:weq))-([6](#fem:sys:wT:ex:Teq)) has\n",
    "only a one-way coupling: $T$ depends on $w$, but $w$ does not depend on\n",
    "$T$. Hence, we can solve ([5](#fem:sys:wT:ex:weq)) with respect\n",
    "to $w$ and then ([6](#fem:sys:wT:ex:Teq)) with respect to $T$.\n",
    "Some may argue that this is not a real system of PDEs, but just two scalar\n",
    "PDEs. Nevertheless, the one-way coupling\n",
    "is convenient when comparing different variational forms\n",
    "and different implementations.\n",
    "\n",
    "# Identical function spaces for the unknowns\n",
    "\n",
    "Let us first apply the same function space $V$ for $w$ and $T$\n",
    "(or more precisely, $w\\in V$ and $T-T_0 \\in V$).\n",
    "With"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "V = \\hbox{span}\\{{\\psi}_0(x,y),\\ldots,{\\psi}_N(x,y)\\},\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "we write"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:sys:wT:ex:sum\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "w = \\sum_{j=0}^N c^{(w)}_j {\\psi}_j,\\quad T = T_0 + \\sum_{j=0}^N c^{(T)}_j\n",
    "{\\psi}_j{\\thinspace .}\n",
    "\\label{fem:sys:wT:ex:sum} \\tag{7}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that $w$ and $T$ in ([5](#fem:sys:wT:ex:weq))-([6](#fem:sys:wT:ex:Teq))\n",
    "denote the exact solution of the PDEs, while $w$ and $T$\n",
    "in ([7](#fem:sys:wT:ex:sum)) are the discrete functions that approximate\n",
    "the exact solution. It should be clear from the context whether a\n",
    "symbol means the exact or approximate solution, but when we need both\n",
    "at the same time, we use a subscript e to denote the exact solution.\n",
    "\n",
    "## Variational form of each individual PDE\n",
    "\n",
    "Inserting the expansions ([7](#fem:sys:wT:ex:sum))\n",
    "in the governing PDEs, results in a residual in each equation,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:sys:wT:ex:weq:R\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "R_w = \\mu \\nabla^2 w + \\beta,\n",
    "\\label{fem:sys:wT:ex:weq:R} \\tag{8}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:sys:wT:ex:Teq:R\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "R_T = \\kappa\\nabla^2 T + \\mu ||\\nabla w||^2\n",
    "{\\thinspace .}\n",
    "\\label{fem:sys:wT:ex:Teq:R} \\tag{9}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A Galerkin method demands $R_w$ and $R_T$ do be orthogonal to $V$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "\\int_\\Omega R_w v {\\, \\mathrm{d}x} &=0\\quad\\forall v\\in V,\\\\\n",
    "\\int_\\Omega R_T v {\\, \\mathrm{d}x} &=0\\quad\\forall v\\in V\n",
    "{\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because of the Dirichlet conditions, $v=0$ on $\\partial\\Omega$.\n",
    "We integrate the Laplace terms by parts and note that the boundary terms\n",
    "vanish since $v=0$ on $\\partial\\Omega$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:sys:wT:ex:w:vf1\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\int_\\Omega \\mu \\nabla w\\cdot\\nabla v {\\, \\mathrm{d}x} = \\int_\\Omega \\beta v{\\, \\mathrm{d}x}\n",
    "\\quad\\forall v\\in V,\n",
    "\\label{fem:sys:wT:ex:w:vf1} \\tag{10}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:sys:wT:ex:T:vf1\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "\\int_\\Omega \\kappa \\nabla T\\cdot\\nabla v {\\, \\mathrm{d}x} = \\int_\\Omega \\mu\n",
    "\\nabla w\\cdot\\nabla w\\, v{\\, \\mathrm{d}x} \\quad\\forall v\\in V\n",
    "\\label{fem:sys:wT:ex:T:vf1} \\tag{11}\n",
    "{\\thinspace .}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The equation $R_w$ in ([8](#fem:sys:wT:ex:weq:R)) is linear\n",
    "in $w$, while the equation $R_T$ in ([9](#fem:sys:wT:ex:Teq:R))\n",
    "is linear in $T$ and nonlinear in $w$.\n",
    "\n",
    "## Compound scalar variational form\n",
    "\n",
    "The alternative way of deriving the variational from is to\n",
    "introduce a test vector function $\\boldsymbol{v}\\in\\boldsymbol{V} = V\\times V$ and take\n",
    "the inner product of $\\boldsymbol{v}$ and the residuals, integrated over the domain:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\int_{\\Omega} (R_w, R_T)\\cdot\\boldsymbol{v} {\\, \\mathrm{d}x} = 0\\quad\\forall\\boldsymbol{v}\\in\\boldsymbol{V}\n",
    "{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With $\\boldsymbol{v} = (v_0,v_1)$ we get"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\int_{\\Omega} (R_w v_0 + R_T v_1) {\\, \\mathrm{d}x} = 0\\quad\\forall\\boldsymbol{v}\\in\\boldsymbol{V}\n",
    "{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Integrating the Laplace terms by parts results in"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:sys:wT:ex:wT:vf2\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\int_\\Omega (\\mu\\nabla w\\cdot\\nabla v_0 + \\kappa\\nabla T\\cdot\\nabla v_1){\\, \\mathrm{d}x}\n",
    "= \\int_\\Omega (\\beta v_0 + \\mu\\nabla w\\cdot\\nabla w\\, v_1){\\, \\mathrm{d}x},\n",
    "\\quad\\forall \\boldsymbol{v}\\in\\boldsymbol{V}\n",
    "{\\thinspace .}\n",
    "\\label{fem:sys:wT:ex:wT:vf2} \\tag{12}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Choosing $v_0=v$ and $v_1=0$ gives the variational form\n",
    "([10](#fem:sys:wT:ex:w:vf1)), while $v_0=0$ and $v_1=v$ gives\n",
    "([11](#fem:sys:wT:ex:T:vf1)).\n",
    "\n",
    "With the inner product notation, $(p,q) = \\int_\\Omega pq{\\, \\mathrm{d}x}$, we\n",
    "can alternatively write ([10](#fem:sys:wT:ex:w:vf1)) and\n",
    "([11](#fem:sys:wT:ex:T:vf1)) as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    " (\\mu\\nabla w,\\nabla v) &= (\\beta, v)\n",
    "\\quad\\forall v\\in V,\\\\\n",
    "(\\kappa \\nabla T,\\nabla v) &= (\\mu\\nabla w\\cdot\\nabla w, v)\\quad\\forall v\\in V,\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or since $\\mu$ and $\\kappa$ are considered constant,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:sys:wT:ex:w:vf1i\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\mu (\\nabla w,\\nabla v) = (\\beta, v)\n",
    "\\quad\\forall v\\in V,\n",
    "\\label{fem:sys:wT:ex:w:vf1i} \\tag{13}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:sys:wT:ex:T:vf1i\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "\\kappa(\\nabla T,\\nabla v) = \\mu(\\nabla w\\cdot\\nabla w, v)\\quad\\forall v\\in V\n",
    "\\label{fem:sys:wT:ex:T:vf1i} \\tag{14}\n",
    "{\\thinspace .}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the left-hand side of ([13](#fem:sys:wT:ex:w:vf1i)) is\n",
    "again linear in $w$, the left-hand side\n",
    "of ([14](#fem:sys:wT:ex:T:vf1i)) is linear in $T$\n",
    "and the nonlinearity of $w$ appears in the right-hand side\n",
    "of  ([14](#fem:sys:wT:ex:T:vf1i))\n",
    "\n",
    "\n",
    "## Decoupled linear systems\n",
    "\n",
    "The linear systems governing the coefficients $c_j^{(w)}$ and\n",
    "$c_j^{(T)}$, $j=0,\\ldots,N$, are derived by inserting the\n",
    "expansions ([7](#fem:sys:wT:ex:sum)) in ([10](#fem:sys:wT:ex:w:vf1))\n",
    "and ([11](#fem:sys:wT:ex:T:vf1)), and choosing $v={\\psi}_i$ for\n",
    "$i=0,\\ldots,N$. The result becomes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:sys:wT:ex:linsys:w1\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\sum_{j=0}^N A^{(w)}_{i,j} c^{(w)}_j = b_i^{(w)},\\quad i=0,\\ldots,N,\n",
    "\\label{fem:sys:wT:ex:linsys:w1} \\tag{15}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:sys:wT:ex:linsys:T1\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "\\sum_{j=0}^N A^{(T)}_{i,j} c^{(T)}_j = b_i^{(T)},\\quad i=0,\\ldots,N,\n",
    "\\label{fem:sys:wT:ex:linsys:T1} \\tag{16}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto2\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "A^{(w)}_{i,j} = \\mu(\\nabla {\\psi}_j,\\nabla {\\psi}_i),\n",
    "\\label{_auto2} \\tag{17}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto3\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "b_i^{(w)} = (\\beta, {\\psi}_i),\n",
    "\\label{_auto3} \\tag{18}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto4\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "A^{(T)}_{i,j} = \\kappa(\\nabla {\\psi}_j,\\nabla {\\psi}_i),\n",
    "\\label{_auto4} \\tag{19}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto5\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "b_i^{(T)} = \\mu((\\sum_j c^{(w)}_j\\nabla{\\psi}_j)\\cdot (\\sum_k\n",
    "c^{(w)}_k\\nabla{\\psi}_k), {\\psi}_i)\n",
    "{\\thinspace .}\n",
    "\\label{_auto5} \\tag{20}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It can also be instructive to write the linear systems using matrices\n",
    "and vectors. Define $K$ as the matrix corresponding to the Laplace\n",
    "operator $\\nabla^2$. That is, $K_{i,j} = (\\nabla {\\psi}_j,\\nabla {\\psi}_i)$.\n",
    "Let us introduce the vectors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "b^{(w)} &= (b_0^{(w)},\\ldots,b_{N}^{(w)}),\\\\\n",
    "b^{(T)} &= (b_0^{(T)},\\ldots,b_{N}^{(T)}),\\\\\n",
    "c^{(w)} &= (c_0^{(w)},\\ldots,c_{N}^{(w)}),\\\\\n",
    "c^{(T)} &= (c_0^{(T)},\\ldots,c_{N}^{(T)}){\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The system ([15](#fem:sys:wT:ex:linsys:w1))-([16](#fem:sys:wT:ex:linsys:T1))\n",
    "can now be expressed in matrix-vector form as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto6\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\mu K c^{(w)} = b^{(w)},\n",
    "\\label{_auto6} \\tag{21}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto7\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "\\kappa K c^{(T)} = b^{(T)}{\\thinspace .}\n",
    "\\label{_auto7} \\tag{22}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can solve the first system for $c^{(w)}$, and then\n",
    "the right-hand side $b^{(T)}$ is known such that we can\n",
    "solve the second system for $c^{(T)}$. Hence, the\n",
    "decoupling of the unknowns $w$ and $T$ reduces the\n",
    "system of nonlinear PDEs to two linear PDEs.\n",
    "\n",
    "\n",
    "## Coupled linear systems\n",
    "\n",
    "Despite the fact that $w$ can be computed first, without knowing $T$,\n",
    "we shall now pretend that $w$ and $T$ enter a two-way coupling such\n",
    "that we need to derive the\n",
    "algebraic equations as *one system* for all the unknowns\n",
    "$c_j^{(w)}$ and $c_j^{(T)}$, $j=0,\\ldots,N$. This system is\n",
    "nonlinear in $c_j^{(w)}$ because of the $\\nabla w\\cdot\\nabla w$ product.\n",
    "To remove this nonlinearity, imagine that we introduce an iteration\n",
    "method where we replace $\\nabla w\\cdot\\nabla w$ by\n",
    "$\\nabla w_{-}\\cdot\\nabla w$, $w_{-}$ being the $w$\n",
    "computed in the previous iteration. Then the term\n",
    "$\\nabla w_{-}\\cdot\\nabla w$ is linear in $w$ since $w_{-}$ is\n",
    "known. The total linear system becomes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:sys:wT:ex:linsys:w2\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\sum_{j=0}^N A^{(w,w)}_{i,j} c^{(w)}_j + \\sum_{j=0}^N A^{(w,T)}_{i,j} c^{(T)}_j\n",
    "= b_i^{(w)},\\quad i=0,\\ldots,N,\n",
    "\\label{fem:sys:wT:ex:linsys:w2} \\tag{23}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:sys:wT:ex:linsys:T2\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "\\sum_{j=0}^N A^{(T,w)}_{i,j} c^{(w)}_j + \\sum_{j=0}^N A^{(T,T)}_{i,j} c^{(T)}_j = b_i^{(T)},\\quad i=0,\\ldots,N,\n",
    "\\label{fem:sys:wT:ex:linsys:T2} \\tag{24}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto8\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "A^{(w,w)}_{i,j} = \\mu(\\nabla {\\psi}_j,\\nabla {\\psi}_i),\n",
    "\\label{_auto8} \\tag{25}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto9\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "A^{(w,T)}_{i,j} = 0,\n",
    "\\label{_auto9} \\tag{26}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto10\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "b_i^{(w)} = (\\beta, {\\psi}_i),\n",
    "\\label{_auto10} \\tag{27}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto11\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "A^{(w,T)}_{i,j} = \\mu((\\nabla w_{-})\\cdot\\nabla{\\psi}_j), {\\psi}_i),\n",
    "\\label{_auto11} \\tag{28}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto12\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "A^{(T,T)}_{i,j} = \\kappa(\\nabla {\\psi}_j,\\nabla {\\psi}_i),\n",
    "\\label{_auto12} \\tag{29}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto13\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "b_i^{(T)} = 0\n",
    "{\\thinspace .}\n",
    "\\label{_auto13} \\tag{30}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This system can alternatively be written in matrix-vector form as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto14\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\mu K c^{(w)} = b^{(w)},\n",
    "\\label{_auto14} \\tag{31}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto15\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "L c^{(w)} + \\kappa K c^{(T)}  =0,\n",
    "\\label{_auto15} \\tag{32}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "with $L$ as the matrix from the $\\nabla w_{-}\\cdot\\nabla$ operator:\n",
    "$L_{i,j} = A^{(w,T)}_{i,j}$. The matrix $K$ is $K_{i,j} =\n",
    "A^{(w,w)}_{i,j} = A^{(T,T)}_{i,j}$.\n",
    "\n",
    "The matrix-vector equations are often conveniently written in block form:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\left(\\begin{array}{cc}\n",
    "\\mu K & 0\\\\\n",
    "L & \\kappa K\n",
    "\\end{array}\\right)\n",
    "\\left(\\begin{array}{c}\n",
    "c^{(w)}\\\\\n",
    "c^{(T)}\n",
    "\\end{array}\\right) =\n",
    "\\left(\\begin{array}{c}\n",
    "b^{(w)}\\\\\n",
    "0\n",
    "\\end{array}\\right),\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that in the general case where all unknowns enter all equations,\n",
    "we have to solve the compound system\n",
    "([23](#fem:sys:wT:ex:linsys:w2))-([24](#fem:sys:wT:ex:linsys:T2)) since\n",
    "then we cannot utilize the special property that\n",
    "([15](#fem:sys:wT:ex:linsys:w1)) does not involve $T$ and can be solved\n",
    "first.\n",
    "\n",
    "When the viscosity depends on the temperature, the\n",
    "$\\mu\\nabla^2w$ term must be replaced by $\\nabla\\cdot (\\mu(T)\\nabla w)$,\n",
    "and then $T$ enters the equation for $w$. Now we have a two-way coupling\n",
    "since both equations contain $w$ and $T$ and therefore\n",
    "must be solved simultaneously.\n",
    "The equation $\\nabla\\cdot (\\mu(T)\\nabla w)=-\\beta$ is nonlinear,\n",
    "and if some iteration procedure is invoked, where we use a previously\n",
    "computed $T_{-}$ in the viscosity ($\\mu(T_{-})$), the coefficient is known,\n",
    "and the equation involves only one unknown, $w$. In that case we are\n",
    "back to the one-way coupled set of PDEs.\n",
    "\n",
    "\n",
    "We may also formulate our PDE system as a vector equation. To this end,\n",
    "we introduce the vector of unknowns $\\boldsymbol{u} = (u^{(0)},u^{(1)})$,\n",
    "where $u^{(0)}=w$ and $u^{(1)}=T$. We then have"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\nabla^2 \\boldsymbol{u} = \\left(\\begin{array}{cc}\n",
    "-{\\mu}^{-1}{\\beta}\\\\\n",
    "-{\\kappa}^{-1}\\mu \\nabla u^{(0)}\\cdot\\nabla u^{(0)}\n",
    "\\end{array}\\right){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Different function spaces for the unknowns\n",
    "\n",
    "\n",
    "It is easy to generalize the previous formulation to the case where\n",
    "$w\\in V^{(w)}$ and $T\\in V^{(T)}$, where $V^{(w)}$ and $V^{(T)}$\n",
    "can be different spaces with different numbers of degrees of freedom.\n",
    "For example, we may use quadratic basis functions for $w$ and linear\n",
    "for $T$. Approximation of the unknowns by different finite element\n",
    "spaces is known as *mixed finite element methods*.\n",
    "\n",
    "We write"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "V^{(w)} &= \\hbox{span}\\{{\\psi}_0^{(w)},\\ldots,{\\psi}_{N_w}^{(w)}\\},\\\\\n",
    "V^{(T)} &= \\hbox{span}\\{{\\psi}_0^{(T)},\\ldots,{\\psi}_{N_T}^{(T)}\\}\n",
    "{\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to\n",
    "multiply ([5](#fem:sys:wT:ex:weq)) by a test function $v^{(w)}\\in V^{(w)}$\n",
    "and ([6](#fem:sys:wT:ex:Teq)) by a $v^{(T)}\\in V^{(T)}$, integrate by\n",
    "parts and arrive at"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:sys:wT:ex:w:vf3\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\int_\\Omega \\mu \\nabla w\\cdot\\nabla v^{(w)} {\\, \\mathrm{d}x} = \\int_\\Omega \\beta v^{(w)}{\\, \\mathrm{d}x}\n",
    "\\quad\\forall v^{(w)}\\in V^{(w)},\n",
    "\\label{fem:sys:wT:ex:w:vf3} \\tag{33}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:sys:wT:ex:T:vf3\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "\\int_\\Omega \\kappa \\nabla T\\cdot\\nabla v^{(T)} {\\, \\mathrm{d}x} = \\int_\\Omega \\mu\n",
    "\\nabla w\\cdot\\nabla w\\, v^{(T)}{\\, \\mathrm{d}x} \\quad\\forall v^{(T)}\\in V^{(T)}\n",
    "\\label{fem:sys:wT:ex:T:vf3} \\tag{34}\n",
    "{\\thinspace .}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The compound scalar variational formulation applies a test vector function\n",
    "$\\boldsymbol{v} = (v^{(w)}, v^{(T)})$ and reads"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:sys:wT:ex:wT:vf3\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\int_\\Omega (\\mu\\nabla w\\cdot\\nabla v^{(w)} +\n",
    "\\kappa\\nabla T\\cdot\\nabla v^{(T)}){\\, \\mathrm{d}x}\n",
    "= \\int_\\Omega (\\beta v^{(w)} + \\mu\\nabla w\\cdot\\nabla w\\, v^{(T)}){\\, \\mathrm{d}x},\n",
    "\\label{fem:sys:wT:ex:wT:vf3} \\tag{35}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "valid $\\forall \\boldsymbol{v}\\in\\boldsymbol{V} = V^{(w)}\\times V^{(T)}$.\n",
    "\n",
    "As earlier, we may decoupled the system in terms\n",
    "of two linear PDEs as we did with\n",
    "([15](#fem:sys:wT:ex:linsys:w1))-([16](#fem:sys:wT:ex:linsys:T1))\n",
    "or linearize the coupled system by introducing the previous\n",
    "iterate $w_{-}$ as in\n",
    "([23](#fem:sys:wT:ex:linsys:w2))-([24](#fem:sys:wT:ex:linsys:T2)).\n",
    "However, we need to distinguish between ${\\psi}_i^{(w)}$\n",
    "and ${\\psi}_i^{(T)}$, and the range in the sums over $j$\n",
    "must match the number of degrees of freedom in the spaces $V^{(w)}$\n",
    "and $V^{(T)}$. The formulas become"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:sys:wT:ex:linsys:w1:mixed\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\sum_{j=0}^{N_w} A^{(w)}_{i,j} c^{(w)}_j = b_i^{(w)},\\quad i=0,\\ldots,N_w,\n",
    "\\label{fem:sys:wT:ex:linsys:w1:mixed} \\tag{36}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:sys:wT:ex:linsys:T1:mixed\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "\\sum_{j=0}^{N_T} A^{(T)}_{i,j} c^{(T)}_j = b_i^{(T)},\\quad i=0,\\ldots,N_T,\n",
    "\\label{fem:sys:wT:ex:linsys:T1:mixed} \\tag{37}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto16\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "A^{(w)}_{i,j} = \\mu(\\nabla {\\psi}_j^{(w)},\\nabla {\\psi}_i^{(w)}),\n",
    "\\label{_auto16} \\tag{38}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto17\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "b_i^{(w)} = (\\beta, {\\psi}_i^{(w)}),\n",
    "\\label{_auto17} \\tag{39}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto18\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "A^{(T)}_{i,j} = \\kappa(\\nabla {\\psi}_j^{(T)},\\nabla {\\psi}_i^{(T)}),\n",
    "\\label{_auto18} \\tag{40}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto19\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "b_i^{(T)} = \\mu(\\sum_{j=0}^{N_w} c^{(w)}_j\\nabla{\\psi}_j^{(w)})\\cdot (\\sum_{k=0}^{N_w}\n",
    "c^{(w)}_k\\nabla{\\psi}_k^{(w)}) , {\\psi}_i^{(T)})\n",
    "{\\thinspace .}\n",
    "\\label{_auto19} \\tag{41}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the case we formulate one compound linear system involving\n",
    "both $c^{(w)}_j$, $j=0,\\ldots,N_w$, and $c^{(T)}_j$, $j=0,\\ldots,N_T$,\n",
    "([23](#fem:sys:wT:ex:linsys:w2))-([24](#fem:sys:wT:ex:linsys:T2))\n",
    "becomes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:sys:wT:ex:linsys:w2b\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\sum_{j=0}^{N_w} A^{(w,w)}_{i,j} c^{(w)}_j + \\sum_{j=0}^{N_T} A^{(w,T)}_{i,j} c^{(T)}_j\n",
    "= b_i^{(w)},\\quad i=0,\\ldots,N_w,\n",
    "\\label{fem:sys:wT:ex:linsys:w2b} \\tag{42}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:sys:wT:ex:linsys:T2b\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "\\sum_{j=0}^{N_w} A^{(T,w)}_{i,j} c^{(w)}_j + \\sum_{j=0}^{N_T} A^{(T,T)}_{i,j} c^{(T)}_j = b_i^{(T)},\\quad i=0,\\ldots,N_T,\n",
    "\\label{fem:sys:wT:ex:linsys:T2b} \\tag{43}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto20\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "A^{(w,w)}_{i,j} = \\mu(\\nabla {\\psi}_j^{(w)},\\nabla {\\psi}_i^{(w)}),\n",
    "\\label{_auto20} \\tag{44}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto21\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "A^{(w,T)}_{i,j} = 0,\n",
    "\\label{_auto21} \\tag{45}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto22\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "b_i^{(w)} = (\\beta, {\\psi}_i^{(w)}),\n",
    "\\label{_auto22} \\tag{46}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto23\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "A^{(w,T)}_{i,j} = \\mu (\\nabla w_{-}\\cdot\\nabla{\\psi}_j^{(w)}), {\\psi}_i^{(T)}),\n",
    "\\label{_auto23} \\tag{47}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto24\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "A^{(T,T)}_{i,j} = \\kappa(\\nabla {\\psi}_j^{(T)},\\nabla {\\psi}_i^{(T)}),\n",
    "\\label{_auto24} \\tag{48}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto25\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "b_i^{(T)} = 0\n",
    "{\\thinspace .}\n",
    "\\label{_auto25} \\tag{49}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we have again performed a linearization by employing a previous iterate $w_{-}$.\n",
    "The corresponding block form"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\left(\\begin{array}{cc}\n",
    "\\mu K^{(w)} & 0\\\\\n",
    "L & \\kappa K^{(T)}\n",
    "\\end{array}\\right)\n",
    "\\left(\\begin{array}{c}\n",
    "c^{(w)}\\\\\n",
    "c^{(T)}\n",
    "\\end{array}\\right) =\n",
    "\\left(\\begin{array}{c}\n",
    "b^{(w)}\\\\\n",
    "0\n",
    "\\end{array}\\right),\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "has square and rectangular block matrices: $K^{(w)}$ is $N_w\\times N_w$,\n",
    "$K^{(T)}$ is $N_T\\times N_T$, while $L$ is $N_T\\times N_w$,\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}