Variational forms for systems of PDEs¶

Many mathematical models involve $m+1$ unknown functions governed by a system of $m+1$ differential equations. In abstract form we may denote the unknowns by $u^{(0)},\ldots, u^{(m)}$ and write the governing equations as

$$\begin{align*} \mathcal{L}_0(u^{(0)},\ldots,u^{(m)}) &= 0,\\ &\vdots\\ \mathcal{L}_{m}(u^{(0)},\ldots,u^{(m)}) &= 0, \end{align*}$$

where $\mathcal{L}_i$ is some differential operator defining differential equation number $i$.

Variational forms¶

There are basically two ways of formulating a variational form for a system of differential equations. The first method treats each equation independently as a scalar equation, while the other method views the total system as a vector equation with a vector function as unknown.

Sequence of scalar PDEs formulation¶

Let us start with the approach that treats one equation at a time. We multiply equation number $i$ by some test function $v^{(i)}\in V^{(i)}$ and integrate over the domain:

$$\begin{equation} \int_\Omega \mathcal{L}^{(0)}(u^{(0)},\ldots,u^{(m)}) v^{(0)}{\, \mathrm{d}x} = 0, \label{fem:sys:vform:1by1a} \tag{1} \end{equation}$$

$$\begin{equation} \vdots \label{_auto1} \tag{2} \end{equation}$$

$$\begin{equation} \int_\Omega \mathcal{L}^{(m)}(u^{(0)},\ldots,u^{(m)}) v^{(m)}{\, \mathrm{d}x} = 0 \label{fem:sys:vform:1by1b} \tag{3} {\thinspace .} \end{equation}$$

Terms with second-order derivatives may be integrated by parts, with Neumann conditions inserted in boundary integrals. Let

$$V^{(i)} = \hbox{span}\{{\psi}_0^{(i)},\ldots,{\psi}_{N_i}^{(i)}\},$$

such that

$$u^{(i)} = B^{(i)}(\boldsymbol{x}) + \sum_{j=0}^{N_i} c_j^{(i)} {\psi}_j^{(i)}(\boldsymbol{x}),$$

where $B^{(i)}$ is a boundary function to handle nonzero Dirichlet conditions. Observe that different unknowns may live in different spaces with different basis functions and numbers of degrees of freedom.

From the $m$ equations in the variational forms we can derive $m$ coupled systems of algebraic equations for the $\Pi_{i=0}^{m} N_i$ unknown coefficients $c_j^{(i)}$, $j=0,\ldots,N_i$, $i=0,\ldots,m$.

Vector PDE formulation¶

The alternative method for deriving a variational form for a system of differential equations introduces a vector of unknown functions

$$\boldsymbol{u} = (u^{(0)},\ldots,u^{(m)}),$$

a vector of test functions

$$\boldsymbol{v} = (v^{(0)},\ldots,v^{(m)}),$$

with

$$\boldsymbol{u}, \boldsymbol{v} \in \boldsymbol{V} = V^{(0)}\times \cdots \times V^{(m)} {\thinspace .}$$

With nonzero Dirichlet conditions, we have a vector $\boldsymbol{B} = (B^{(0)},\ldots,B^{(m)})$ with boundary functions and then it is $\boldsymbol{u} - \boldsymbol{B}$ that lies in $\boldsymbol{V}$, not $\boldsymbol{u}$ itself.

The governing system of differential equations is written

$$\boldsymbol{\mathcal{L}}(\boldsymbol{u} ) = 0,$$

where

$$\boldsymbol{\mathcal{L}}(\boldsymbol{u} ) = (\mathcal{L}^{(0)}(\boldsymbol{u}),\ldots, \mathcal{L}^{(m)}(\boldsymbol{u})) {\thinspace .}$$

The variational form is derived by taking the inner product of the vector of equations and the test function vector:

$$\begin{equation} \int_\Omega \boldsymbol{\mathcal{L}}(\boldsymbol{u} )\cdot\boldsymbol{v} = 0\quad\forall\boldsymbol{v}\in\boldsymbol{V}{\thinspace .} \label{fem:sys:vform:inner} \tag{4} \end{equation}$$

Observe that (4) is one scalar equation. To derive systems of algebraic equations for the unknown coefficients in the expansions of the unknown functions, one chooses $m$ linearly independent $\boldsymbol{v}$ vectors to generate $m$ independent variational forms from (4). The particular choice $\boldsymbol{v} = (v^{(0)},0,\ldots,0)$ recovers (1), $\boldsymbol{v} = (0,\ldots,0,v^{(m)})$ recovers (3), and $\boldsymbol{v} = (0,\ldots,0,v^{(i)},0,\ldots,0)$ recovers the variational form number $i$, $\int_\Omega \mathcal{L}^{(i)} v^{(i)}{\, \mathrm{d}x} =0$, in (1)-(3).

A worked example¶

We now consider a specific system of two partial differential equations in two space dimensions:

$$\begin{equation} \mu \nabla^2 w = -\beta, \label{fem:sys:wT:ex:weq} \tag{5} \end{equation}$$

$$\begin{equation} \kappa\nabla^2 T = - \mu ||\nabla w||^2 {\thinspace .} \label{fem:sys:wT:ex:Teq} \tag{6} \end{equation}$$

The unknown functions $w(x,y)$ and $T(x,y)$ are defined in a domain $\Omega$, while $\mu$, $\beta$, and $\kappa$ are given constants. The norm in (6) is the standard Euclidean norm:

$$||\nabla w||^2 = \nabla w\cdot\nabla w = w_x^2 + w_y^2 {\thinspace .}$$

The boundary conditions associated with (5)-(6) are $w=0$ on $\partial\Omega$ and $T=T_0$ on $\partial\Omega$. Each of the equations (5) and (6) needs one condition at each point on the boundary.

The system (5)-(6) arises from fluid flow in a straight pipe, with the $z$ axis in the direction of the pipe. The domain $\Omega$ is a cross section of the pipe, $w$ is the velocity in the $z$ direction, $\mu$ is the viscosity of the fluid, $\beta$ is the pressure gradient along the pipe, $T$ is the temperature, and $\kappa$ is the heat conduction coefficient of the fluid. The equation (5) comes from the Navier-Stokes equations, and (6) follows from the energy equation. The term $- \mu ||\nabla w||^2$ models heating of the fluid due to internal friction.

Observe that the system (5)-(6) has only a one-way coupling: $T$ depends on $w$, but $w$ does not depend on $T$. Hence, we can solve (5) with respect to $w$ and then (6) with respect to $T$. Some may argue that this is not a real system of PDEs, but just two scalar PDEs. Nevertheless, the one-way coupling is convenient when comparing different variational forms and different implementations.

Identical function spaces for the unknowns¶

Let us first apply the same function space $V$ for $w$ and $T$ (or more precisely, $w\in V$ and $T-T_0 \in V$). With

$$V = \hbox{span}\{{\psi}_0(x,y),\ldots,{\psi}_N(x,y)\},$$

we write

$$\begin{equation} w = \sum_{j=0}^N c^{(w)}_j {\psi}_j,\quad T = T_0 + \sum_{j=0}^N c^{(T)}_j {\psi}_j{\thinspace .} \label{fem:sys:wT:ex:sum} \tag{7} \end{equation}$$

Note that $w$ and $T$ in (5)-(6) denote the exact solution of the PDEs, while $w$ and $T$ in (7) are the discrete functions that approximate the exact solution. It should be clear from the context whether a symbol means the exact or approximate solution, but when we need both at the same time, we use a subscript e to denote the exact solution.

Variational form of each individual PDE¶

Inserting the expansions (7) in the governing PDEs, results in a residual in each equation,

$$\begin{equation} R_w = \mu \nabla^2 w + \beta, \label{fem:sys:wT:ex:weq:R} \tag{8} \end{equation}$$

$$\begin{equation} R_T = \kappa\nabla^2 T + \mu ||\nabla w||^2 {\thinspace .} \label{fem:sys:wT:ex:Teq:R} \tag{9} \end{equation}$$

A Galerkin method demands $R_w$ and $R_T$ do be orthogonal to $V$:

$$\begin{align*} \int_\Omega R_w v {\, \mathrm{d}x} &=0\quad\forall v\in V,\\ \int_\Omega R_T v {\, \mathrm{d}x} &=0\quad\forall v\in V {\thinspace .} \end{align*}$$

Because of the Dirichlet conditions, $v=0$ on $\partial\Omega$. We integrate the Laplace terms by parts and note that the boundary terms vanish since $v=0$ on $\partial\Omega$:

$$\begin{equation} \int_\Omega \mu \nabla w\cdot\nabla v {\, \mathrm{d}x} = \int_\Omega \beta v{\, \mathrm{d}x} \quad\forall v\in V, \label{fem:sys:wT:ex:w:vf1} \tag{10} \end{equation}$$

$$\begin{equation} \int_\Omega \kappa \nabla T\cdot\nabla v {\, \mathrm{d}x} = \int_\Omega \mu \nabla w\cdot\nabla w\, v{\, \mathrm{d}x} \quad\forall v\in V \label{fem:sys:wT:ex:T:vf1} \tag{11} {\thinspace .} \end{equation}$$

The equation $R_w$ in (8) is linear in $w$, while the equation $R_T$ in (9) is linear in $T$ and nonlinear in $w$.

Compound scalar variational form¶

The alternative way of deriving the variational from is to introduce a test vector function $\boldsymbol{v}\in\boldsymbol{V} = V\times V$ and take the inner product of $\boldsymbol{v}$ and the residuals, integrated over the domain:

$$\int_{\Omega} (R_w, R_T)\cdot\boldsymbol{v} {\, \mathrm{d}x} = 0\quad\forall\boldsymbol{v}\in\boldsymbol{V} {\thinspace .}$$

With $\boldsymbol{v} = (v_0,v_1)$ we get

$$\int_{\Omega} (R_w v_0 + R_T v_1) {\, \mathrm{d}x} = 0\quad\forall\boldsymbol{v}\in\boldsymbol{V} {\thinspace .}$$

Integrating the Laplace terms by parts results in

$$\begin{equation} \int_\Omega (\mu\nabla w\cdot\nabla v_0 + \kappa\nabla T\cdot\nabla v_1){\, \mathrm{d}x} = \int_\Omega (\beta v_0 + \mu\nabla w\cdot\nabla w\, v_1){\, \mathrm{d}x}, \quad\forall \boldsymbol{v}\in\boldsymbol{V} {\thinspace .} \label{fem:sys:wT:ex:wT:vf2} \tag{12} \end{equation}$$

Choosing $v_0=v$ and $v_1=0$ gives the variational form (10), while $v_0=0$ and $v_1=v$ gives (11).

With the inner product notation, $(p,q) = \int_\Omega pq{\, \mathrm{d}x}$, we can alternatively write (10) and (11) as

$$\begin{align*} (\mu\nabla w,\nabla v) &= (\beta, v) \quad\forall v\in V,\\ (\kappa \nabla T,\nabla v) &= (\mu\nabla w\cdot\nabla w, v)\quad\forall v\in V, \end{align*}$$

or since $\mu$ and $\kappa$ are considered constant,

$$\begin{equation} \mu (\nabla w,\nabla v) = (\beta, v) \quad\forall v\in V, \label{fem:sys:wT:ex:w:vf1i} \tag{13} \end{equation}$$

$$\begin{equation} \kappa(\nabla T,\nabla v) = \mu(\nabla w\cdot\nabla w, v)\quad\forall v\in V \label{fem:sys:wT:ex:T:vf1i} \tag{14} {\thinspace .} \end{equation}$$

Note that the left-hand side of (13) is again linear in $w$, the left-hand side of (14) is linear in $T$ and the nonlinearity of $w$ appears in the right-hand side of (14)

Decoupled linear systems¶

The linear systems governing the coefficients $c_j^{(w)}$ and $c_j^{(T)}$, $j=0,\ldots,N$, are derived by inserting the expansions (7) in (10) and (11), and choosing $v={\psi}_i$ for $i=0,\ldots,N$. The result becomes

$$\begin{equation} \sum_{j=0}^N A^{(w)}_{i,j} c^{(w)}_j = b_i^{(w)},\quad i=0,\ldots,N, \label{fem:sys:wT:ex:linsys:w1} \tag{15} \end{equation}$$

$$\begin{equation} \sum_{j=0}^N A^{(T)}_{i,j} c^{(T)}_j = b_i^{(T)},\quad i=0,\ldots,N, \label{fem:sys:wT:ex:linsys:T1} \tag{16} \end{equation}$$

$$\begin{equation} A^{(w)}_{i,j} = \mu(\nabla {\psi}_j,\nabla {\psi}_i), \label{_auto2} \tag{17} \end{equation}$$

$$\begin{equation} b_i^{(w)} = (\beta, {\psi}_i), \label{_auto3} \tag{18} \end{equation}$$

$$\begin{equation} A^{(T)}_{i,j} = \kappa(\nabla {\psi}_j,\nabla {\psi}_i), \label{_auto4} \tag{19} \end{equation}$$

$$\begin{equation} b_i^{(T)} = \mu((\sum_j c^{(w)}_j\nabla{\psi}_j)\cdot (\sum_k c^{(w)}_k\nabla{\psi}_k), {\psi}_i) {\thinspace .} \label{_auto5} \tag{20} \end{equation}$$

It can also be instructive to write the linear systems using matrices and vectors. Define $K$ as the matrix corresponding to the Laplace operator $\nabla^2$. That is, $K_{i,j} = (\nabla {\psi}_j,\nabla {\psi}_i)$. Let us introduce the vectors

$$\begin{align*} b^{(w)} &= (b_0^{(w)},\ldots,b_{N}^{(w)}),\\ b^{(T)} &= (b_0^{(T)},\ldots,b_{N}^{(T)}),\\ c^{(w)} &= (c_0^{(w)},\ldots,c_{N}^{(w)}),\\ c^{(T)} &= (c_0^{(T)},\ldots,c_{N}^{(T)}){\thinspace .} \end{align*}$$

The system (15)-(16) can now be expressed in matrix-vector form as

$$\begin{equation} \mu K c^{(w)} = b^{(w)}, \label{_auto6} \tag{21} \end{equation}$$

$$\begin{equation} \kappa K c^{(T)} = b^{(T)}{\thinspace .} \label{_auto7} \tag{22} \end{equation}$$

We can solve the first system for $c^{(w)}$, and then the right-hand side $b^{(T)}$ is known such that we can solve the second system for $c^{(T)}$. Hence, the decoupling of the unknowns $w$ and $T$ reduces the system of nonlinear PDEs to two linear PDEs.

Coupled linear systems¶

Despite the fact that $w$ can be computed first, without knowing $T$, we shall now pretend that $w$ and $T$ enter a two-way coupling such that we need to derive the algebraic equations as one system for all the unknowns $c_j^{(w)}$ and $c_j^{(T)}$, $j=0,\ldots,N$. This system is nonlinear in $c_j^{(w)}$ because of the $\nabla w\cdot\nabla w$ product. To remove this nonlinearity, imagine that we introduce an iteration method where we replace $\nabla w\cdot\nabla w$ by $\nabla w_{-}\cdot\nabla w$, $w_{-}$ being the $w$ computed in the previous iteration. Then the term $\nabla w_{-}\cdot\nabla w$ is linear in $w$ since $w_{-}$ is known. The total linear system becomes

$$\begin{equation} \sum_{j=0}^N A^{(w,w)}_{i,j} c^{(w)}_j + \sum_{j=0}^N A^{(w,T)}_{i,j} c^{(T)}_j = b_i^{(w)},\quad i=0,\ldots,N, \label{fem:sys:wT:ex:linsys:w2} \tag{23} \end{equation}$$

$$\begin{equation} \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, \label{fem:sys:wT:ex:linsys:T2} \tag{24} \end{equation}$$

$$\begin{equation} A^{(w,w)}_{i,j} = \mu(\nabla {\psi}_j,\nabla {\psi}_i), \label{_auto8} \tag{25} \end{equation}$$

$$\begin{equation} A^{(w,T)}_{i,j} = 0, \label{_auto9} \tag{26} \end{equation}$$

$$\begin{equation} b_i^{(w)} = (\beta, {\psi}_i), \label{_auto10} \tag{27} \end{equation}$$

$$\begin{equation} A^{(w,T)}_{i,j} = \mu((\nabla w_{-})\cdot\nabla{\psi}_j), {\psi}_i), \label{_auto11} \tag{28} \end{equation}$$

$$\begin{equation} A^{(T,T)}_{i,j} = \kappa(\nabla {\psi}_j,\nabla {\psi}_i), \label{_auto12} \tag{29} \end{equation}$$

$$\begin{equation} b_i^{(T)} = 0 {\thinspace .} \label{_auto13} \tag{30} \end{equation}$$

This system can alternatively be written in matrix-vector form as

$$\begin{equation} \mu K c^{(w)} = b^{(w)}, \label{_auto14} \tag{31} \end{equation}$$

$$\begin{equation} L c^{(w)} + \kappa K c^{(T)} =0, \label{_auto15} \tag{32} \end{equation}$$

with $L$ as the matrix from the $\nabla w_{-}\cdot\nabla$ operator: $L_{i,j} = A^{(w,T)}_{i,j}$. The matrix $K$ is $K_{i,j} = A^{(w,w)}_{i,j} = A^{(T,T)}_{i,j}$.

The matrix-vector equations are often conveniently written in block form:

$$\left(\begin{array}{cc} \mu K & 0\\ L & \kappa K \end{array}\right) \left(\begin{array}{c} c^{(w)}\\ c^{(T)} \end{array}\right) = \left(\begin{array}{c} b^{(w)}\\ 0 \end{array}\right),$$

Note that in the general case where all unknowns enter all equations, we have to solve the compound system (23)-(24) since then we cannot utilize the special property that (15) does not involve $T$ and can be solved first.

When the viscosity depends on the temperature, the $\mu\nabla^2w$ term must be replaced by $\nabla\cdot (\mu(T)\nabla w)$, and then $T$ enters the equation for $w$. Now we have a two-way coupling since both equations contain $w$ and $T$ and therefore must be solved simultaneously. The equation $\nabla\cdot (\mu(T)\nabla w)=-\beta$ is nonlinear, and if some iteration procedure is invoked, where we use a previously computed $T_{-}$ in the viscosity ($\mu(T_{-})$), the coefficient is known, and the equation involves only one unknown, $w$. In that case we are back to the one-way coupled set of PDEs.

We may also formulate our PDE system as a vector equation. To this end, we introduce the vector of unknowns $\boldsymbol{u} = (u^{(0)},u^{(1)})$, where $u^{(0)}=w$ and $u^{(1)}=T$. We then have

$$\nabla^2 \boldsymbol{u} = \left(\begin{array}{cc} -{\mu}^{-1}{\beta}\\ -{\kappa}^{-1}\mu \nabla u^{(0)}\cdot\nabla u^{(0)} \end{array}\right){\thinspace .}$$

Different function spaces for the unknowns¶

It is easy to generalize the previous formulation to the case where $w\in V^{(w)}$ and $T\in V^{(T)}$, where $V^{(w)}$ and $V^{(T)}$ can be different spaces with different numbers of degrees of freedom. For example, we may use quadratic basis functions for $w$ and linear for $T$. Approximation of the unknowns by different finite element spaces is known as mixed finite element methods.

We write

$$\begin{align*} V^{(w)} &= \hbox{span}\{{\psi}_0^{(w)},\ldots,{\psi}_{N_w}^{(w)}\},\\ V^{(T)} &= \hbox{span}\{{\psi}_0^{(T)},\ldots,{\psi}_{N_T}^{(T)}\} {\thinspace .} \end{align*}$$

The next step is to multiply (5) by a test function $v^{(w)}\in V^{(w)}$ and (6) by a $v^{(T)}\in V^{(T)}$, integrate by parts and arrive at

$$\begin{equation} \int_\Omega \mu \nabla w\cdot\nabla v^{(w)} {\, \mathrm{d}x} = \int_\Omega \beta v^{(w)}{\, \mathrm{d}x} \quad\forall v^{(w)}\in V^{(w)}, \label{fem:sys:wT:ex:w:vf3} \tag{33} \end{equation}$$

$$\begin{equation} \int_\Omega \kappa \nabla T\cdot\nabla v^{(T)} {\, \mathrm{d}x} = \int_\Omega \mu \nabla w\cdot\nabla w\, v^{(T)}{\, \mathrm{d}x} \quad\forall v^{(T)}\in V^{(T)} \label{fem:sys:wT:ex:T:vf3} \tag{34} {\thinspace .} \end{equation}$$

The compound scalar variational formulation applies a test vector function $\boldsymbol{v} = (v^{(w)}, v^{(T)})$ and reads

$$\begin{equation} \int_\Omega (\mu\nabla w\cdot\nabla v^{(w)} + \kappa\nabla T\cdot\nabla v^{(T)}){\, \mathrm{d}x} = \int_\Omega (\beta v^{(w)} + \mu\nabla w\cdot\nabla w\, v^{(T)}){\, \mathrm{d}x}, \label{fem:sys:wT:ex:wT:vf3} \tag{35} \end{equation}$$

valid $\forall \boldsymbol{v}\in\boldsymbol{V} = V^{(w)}\times V^{(T)}$.

As earlier, we may decoupled the system in terms of two linear PDEs as we did with (15)-(16) or linearize the coupled system by introducing the previous iterate $w_{-}$ as in (23)-(24). However, we need to distinguish between ${\psi}_i^{(w)}$ and ${\psi}_i^{(T)}$, and the range in the sums over $j$ must match the number of degrees of freedom in the spaces $V^{(w)}$ and $V^{(T)}$. The formulas become

$$\begin{equation} \sum_{j=0}^{N_w} A^{(w)}_{i,j} c^{(w)}_j = b_i^{(w)},\quad i=0,\ldots,N_w, \label{fem:sys:wT:ex:linsys:w1:mixed} \tag{36} \end{equation}$$

$$\begin{equation} \sum_{j=0}^{N_T} A^{(T)}_{i,j} c^{(T)}_j = b_i^{(T)},\quad i=0,\ldots,N_T, \label{fem:sys:wT:ex:linsys:T1:mixed} \tag{37} \end{equation}$$

$$\begin{equation} A^{(w)}_{i,j} = \mu(\nabla {\psi}_j^{(w)},\nabla {\psi}_i^{(w)}), \label{_auto16} \tag{38} \end{equation}$$

$$\begin{equation} b_i^{(w)} = (\beta, {\psi}_i^{(w)}), \label{_auto17} \tag{39} \end{equation}$$

$$\begin{equation} A^{(T)}_{i,j} = \kappa(\nabla {\psi}_j^{(T)},\nabla {\psi}_i^{(T)}), \label{_auto18} \tag{40} \end{equation}$$

$$\begin{equation} b_i^{(T)} = \mu(\sum_{j=0}^{N_w} c^{(w)}_j\nabla{\psi}_j^{(w)})\cdot (\sum_{k=0}^{N_w} c^{(w)}_k\nabla{\psi}_k^{(w)}) , {\psi}_i^{(T)}) {\thinspace .} \label{_auto19} \tag{41} \end{equation}$$

In the case we formulate one compound linear system involving both $c^{(w)}_j$, $j=0,\ldots,N_w$, and $c^{(T)}_j$, $j=0,\ldots,N_T$, (23)-(24) becomes

$$\begin{equation} \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 = b_i^{(w)},\quad i=0,\ldots,N_w, \label{fem:sys:wT:ex:linsys:w2b} \tag{42} \end{equation}$$

$$\begin{equation} \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, \label{fem:sys:wT:ex:linsys:T2b} \tag{43} \end{equation}$$

$$\begin{equation} A^{(w,w)}_{i,j} = \mu(\nabla {\psi}_j^{(w)},\nabla {\psi}_i^{(w)}), \label{_auto20} \tag{44} \end{equation}$$

$$\begin{equation} A^{(w,T)}_{i,j} = 0, \label{_auto21} \tag{45} \end{equation}$$

$$\begin{equation} b_i^{(w)} = (\beta, {\psi}_i^{(w)}), \label{_auto22} \tag{46} \end{equation}$$

$$\begin{equation} A^{(w,T)}_{i,j} = \mu (\nabla w_{-}\cdot\nabla{\psi}_j^{(w)}), {\psi}_i^{(T)}), \label{_auto23} \tag{47} \end{equation}$$

$$\begin{equation} A^{(T,T)}_{i,j} = \kappa(\nabla {\psi}_j^{(T)},\nabla {\psi}_i^{(T)}), \label{_auto24} \tag{48} \end{equation}$$

$$\begin{equation} b_i^{(T)} = 0 {\thinspace .} \label{_auto25} \tag{49} \end{equation}$$

Here, we have again performed a linearization by employing a previous iterate $w_{-}$. The corresponding block form

$$\left(\begin{array}{cc} \mu K^{(w)} & 0\\ L & \kappa K^{(T)} \end{array}\right) \left(\begin{array}{c} c^{(w)}\\ c^{(T)} \end{array}\right) = \left(\begin{array}{c} b^{(w)}\\ 0 \end{array}\right),$$

has square and rectangular block matrices: $K^{(w)}$ is $N_w\times N_w$, $K^{(T)}$ is $N_T\times N_T$, while $L$ is $N_T\times N_w$,