# Variational forms for systems of PDEs
<div id="ch:femsys"></div>

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
<div id="fem:sys:vform"></div>

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:

<!-- Equation labels as ordinary links -->
<div id="fem:sys:vform:1by1a"></div>

$$
\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}
$$

<!-- Equation labels as ordinary links -->
<div id="_auto1"></div>

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

<!-- Equation labels as ordinary links -->
<div id="fem:sys:vform:1by1b"></div>

$$
\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:

<!-- Equation labels as ordinary links -->
<div id="fem:sys:vform:inner"></div>

$$
\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](#fem:sys:vform:inner)) 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](#fem:sys:vform:inner)).  The particular choice $\boldsymbol{v} =
(v^{(0)},0,\ldots,0)$ recovers ([1](#fem:sys:vform:1by1a)), $\boldsymbol{v} =
(0,\ldots,0,v^{(m)})$ recovers ([3](#fem:sys:vform:1by1b)), 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](#fem:sys:vform:1by1a))-([3](#fem:sys:vform:1by1b)).

# A worked example
<div id="fem:sys:uT:ex"></div>

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

<!-- Equation labels as ordinary links -->
<div id="fem:sys:wT:ex:weq"></div>

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

<!-- Equation labels as ordinary links -->
<div id="fem:sys:wT:ex:Teq"></div>

$$
\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](#fem:sys:wT:ex:Teq)) 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](#fem:sys:wT:ex:weq))-([6](#fem:sys:wT:ex:Teq)) are $w=0$ on
$\partial\Omega$ and $T=T_0$ on $\partial\Omega$.
Each of the equations ([5](#fem:sys:wT:ex:weq)) and ([6](#fem:sys:wT:ex:Teq))
needs one condition at each point on the boundary.

The system ([5](#fem:sys:wT:ex:weq))-([6](#fem:sys:wT:ex:Teq)) 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](#fem:sys:wT:ex:weq)) comes from the Navier-Stokes
equations, and ([6](#fem:sys:wT:ex:Teq)) 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](#fem:sys:wT:ex:weq))-([6](#fem:sys:wT:ex:Teq)) has
only a one-way coupling: $T$ depends on $w$, but $w$ does not depend on
$T$. Hence, we can solve ([5](#fem:sys:wT:ex:weq)) with respect
to $w$ and then ([6](#fem:sys:wT:ex:Teq)) 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

<!-- Equation labels as ordinary links -->
<div id="fem:sys:wT:ex:sum"></div>

$$
\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](#fem:sys:wT:ex:weq))-([6](#fem:sys:wT:ex:Teq))
denote the exact solution of the PDEs, while $w$ and $T$
in ([7](#fem:sys:wT:ex:sum)) 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](#fem:sys:wT:ex:sum))
in the governing PDEs, results in a residual in each equation,

<!-- Equation labels as ordinary links -->
<div id="fem:sys:wT:ex:weq:R"></div>

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

<!-- Equation labels as ordinary links -->
<div id="fem:sys:wT:ex:Teq:R"></div>

$$
\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$:

<!-- Equation labels as ordinary links -->
<div id="fem:sys:wT:ex:w:vf1"></div>

$$
\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}
$$

<!-- Equation labels as ordinary links -->
<div id="fem:sys:wT:ex:T:vf1"></div>

$$
\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](#fem:sys:wT:ex:weq:R)) is linear
in $w$, while the equation $R_T$ in ([9](#fem:sys:wT:ex:Teq:R))
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

<!-- Equation labels as ordinary links -->
<div id="fem:sys:wT:ex:wT:vf2"></div>

$$
\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](#fem:sys:wT:ex:w:vf1)), while $v_0=0$ and $v_1=v$ gives
([11](#fem:sys:wT:ex:T:vf1)).

With the inner product notation, $(p,q) = \int_\Omega pq{\, \mathrm{d}x}$, we
can alternatively write ([10](#fem:sys:wT:ex:w:vf1)) and
([11](#fem:sys:wT:ex:T:vf1)) 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,

<!-- Equation labels as ordinary links -->
<div id="fem:sys:wT:ex:w:vf1i"></div>

$$
\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}
$$

<!-- Equation labels as ordinary links -->
<div id="fem:sys:wT:ex:T:vf1i"></div>

$$
\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](#fem:sys:wT:ex:w:vf1i)) is
again linear in $w$, the left-hand side
of ([14](#fem:sys:wT:ex:T:vf1i)) is linear in $T$
and the nonlinearity of $w$ appears in the right-hand side
of  ([14](#fem:sys:wT:ex:T:vf1i))


## 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](#fem:sys:wT:ex:sum)) in ([10](#fem:sys:wT:ex:w:vf1))
and ([11](#fem:sys:wT:ex:T:vf1)), and choosing $v={\psi}_i$ for
$i=0,\ldots,N$. The result becomes

<!-- Equation labels as ordinary links -->
<div id="fem:sys:wT:ex:linsys:w1"></div>

$$
\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}
$$

<!-- Equation labels as ordinary links -->
<div id="fem:sys:wT:ex:linsys:T1"></div>

$$
\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}
$$

<!-- Equation labels as ordinary links -->
<div id="_auto2"></div>

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

<!-- Equation labels as ordinary links -->
<div id="_auto3"></div>

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

<!-- Equation labels as ordinary links -->
<div id="_auto4"></div>

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

<!-- Equation labels as ordinary links -->
<div id="_auto5"></div>

$$
\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](#fem:sys:wT:ex:linsys:w1))-([16](#fem:sys:wT:ex:linsys:T1))
can now be expressed in matrix-vector form as

<!-- Equation labels as ordinary links -->
<div id="_auto6"></div>

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

<!-- Equation labels as ordinary links -->
<div id="_auto7"></div>

$$
\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

<!-- Equation labels as ordinary links -->
<div id="fem:sys:wT:ex:linsys:w2"></div>

$$
\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}
$$

<!-- Equation labels as ordinary links -->
<div id="fem:sys:wT:ex:linsys:T2"></div>

$$
\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}
$$

<!-- Equation labels as ordinary links -->
<div id="_auto8"></div>

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

<!-- Equation labels as ordinary links -->
<div id="_auto9"></div>

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

<!-- Equation labels as ordinary links -->
<div id="_auto10"></div>

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

<!-- Equation labels as ordinary links -->
<div id="_auto11"></div>

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

<!-- Equation labels as ordinary links -->
<div id="_auto12"></div>

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

<!-- Equation labels as ordinary links -->
<div id="_auto13"></div>

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

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

<!-- Equation labels as ordinary links -->
<div id="_auto14"></div>

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

<!-- Equation labels as ordinary links -->
<div id="_auto15"></div>

$$
\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](#fem:sys:wT:ex:linsys:w2))-([24](#fem:sys:wT:ex:linsys:T2)) since
then we cannot utilize the special property that
([15](#fem:sys:wT:ex:linsys:w1)) 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](#fem:sys:wT:ex:weq)) by a test function $v^{(w)}\in V^{(w)}$
and ([6](#fem:sys:wT:ex:Teq)) by a $v^{(T)}\in V^{(T)}$, integrate by
parts and arrive at

<!-- Equation labels as ordinary links -->
<div id="fem:sys:wT:ex:w:vf3"></div>

$$
\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}
$$

<!-- Equation labels as ordinary links -->
<div id="fem:sys:wT:ex:T:vf3"></div>

$$
\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

<!-- Equation labels as ordinary links -->
<div id="fem:sys:wT:ex:wT:vf3"></div>

$$
\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](#fem:sys:wT:ex:linsys:w1))-([16](#fem:sys:wT:ex:linsys:T1))
or linearize the coupled system by introducing the previous
iterate $w_{-}$ as in
([23](#fem:sys:wT:ex:linsys:w2))-([24](#fem:sys:wT:ex:linsys:T2)).
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

<!-- Equation labels as ordinary links -->
<div id="fem:sys:wT:ex:linsys:w1:mixed"></div>

$$
\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}
$$

<!-- Equation labels as ordinary links -->
<div id="fem:sys:wT:ex:linsys:T1:mixed"></div>

$$
\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}
$$

<!-- Equation labels as ordinary links -->
<div id="_auto16"></div>

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

<!-- Equation labels as ordinary links -->
<div id="_auto17"></div>

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

<!-- Equation labels as ordinary links -->
<div id="_auto18"></div>

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

<!-- Equation labels as ordinary links -->
<div id="_auto19"></div>

$$
\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](#fem:sys:wT:ex:linsys:w2))-([24](#fem:sys:wT:ex:linsys:T2))
becomes

<!-- Equation labels as ordinary links -->
<div id="fem:sys:wT:ex:linsys:w2b"></div>

$$
\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}
$$

<!-- Equation labels as ordinary links -->
<div id="fem:sys:wT:ex:linsys:T2b"></div>

$$
\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}
$$

<!-- Equation labels as ordinary links -->
<div id="_auto20"></div>

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

<!-- Equation labels as ordinary links -->
<div id="_auto21"></div>

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

<!-- Equation labels as ordinary links -->
<div id="_auto22"></div>

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

<!-- Equation labels as ordinary links -->
<div id="_auto23"></div>

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

<!-- Equation labels as ordinary links -->
<div id="_auto24"></div>

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

<!-- Equation labels as ordinary links -->
<div id="_auto25"></div>

$$
\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$,



0094abcb2c24b5ef0bedd128418704e408f6b267
ipynb
Jupyter Notebook
7- femsys.ipynb
mbarzegary/finite-element-intro
47ef0a3592b823ae71a874ee35850114f16b6d8b
MIT
2021-01-26T13:18:02.000Z
2022-02-14T15:20:11.000Z
7- femsys.ipynb
mbarzegary/finite-element-intro
47ef0a3592b823ae71a874ee35850114f16b6d8b
MIT
7- femsys.ipynb
mbarzegary/finite-element-intro
47ef0a3592b823ae71a874ee35850114f16b6d8b
MIT
2021-08-05T23:14:15.000Z
2021-10-05T10:22:29.000Z
Qwen/Qwen-72B
1. YES
2. YES