# Systems of differential equations¶

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{split}\mathcal{L}_0(u^{(0)},\ldots,u^{(m)}) &= 0,\\ &\vdots\\ \mathcal{L}_{m}(u^{(0)},\ldots,u^{(m)}) &= 0,\end{split}$

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

## Variational forms (3)¶

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.

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

(1)$\int_\Omega \mathcal{L}^{(0)}(u^{(0)},\ldots,u^{(m)}) v^{(0)}{\, \mathrm{d}x} = 0,$
$\vdots$
(2)$\int_\Omega \mathcal{L}^{(m)}(u^{(0)},\ldots,u^{(m)}) v^{(m)}{\, \mathrm{d}x} = 0$${\thinspace .}$

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

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} = (u^{(0)},\ldots,u^{(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:

(3)$\int_\Omega \boldsymbol{\mathcal{L}}(\boldsymbol{u} )\cdot\boldsymbol{v} = 0\quad\forall\boldsymbol{v}\in\boldsymbol{V}{\thinspace .}$

Observe that (3) 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 (3). The particular choice $$\boldsymbol{v} = (v^{(0)},0,\ldots,0)$$ recovers (1), $$\boldsymbol{v} = (0,\ldots,0,v^{(m)}$$ recovers (2), 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)-(2).

## A worked example¶

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

(4)$\mu \nabla^2 w = -\beta,$
(5)$\kappa\nabla^2 T = - \mu ||\nabla w||^2 {\thinspace .}$

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 (5) is the standard Eucledian norm:

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

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

The system (4)-(5) 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 (4) comes from the Navier-Stokes equations, and (5) follows from the energy equation. The term $$- \mu ||\nabla w||^2$$ models heating of the fluid due to internal friction.

Observe that the system (4)-(5) has only a one-way coupling: $$T$$ depends on $$w$$, but $$w$$ does not depend on $$T$$, because we can solve (4) with respect to $$w$$ and then (5) 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

(6)$w = \sum_{j=0}^N c^{(w)}_j {\psi}_j,\quad T = T_0 + \sum_{j=0}^N c^{(T)}_j {\psi}_j{\thinspace .}$

Note that $$w$$ and $$T$$ in (4)-(5) denote the exact solution of the PDEs, while $$w$$ and $$T$$ (6) 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 (6) in the governing PDEs, results in a residual in each equation,

(7)$R_w = \mu \nabla^2 w + \beta,$
(8)$R_T = \kappa\nabla^2 T + \mu ||\nabla w||^2 {\thinspace .}$

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

$\begin{split}\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{split}$

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

(9)$\int_\Omega \mu \nabla w\cdot\nabla v {\, \mathrm{d}x} = \int_\Omega \beta v{\, \mathrm{d}x} \quad\forall v\in V,$
(10)$\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$${\thinspace .}$

### 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

(11)$\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 .}$

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

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

$\begin{split} (\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{split}$

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

(12)$\mu (\nabla w,\nabla v) = (\beta, v) \quad\forall v\in V,$
(13)$\kappa(\nabla T,\nabla v) = \mu(\nabla w\cdot\nabla w, v)\quad\forall v\in V$${\thinspace .}$

### 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 (6) in (9) and (10), and choosing $$v={\psi}_i$$ for $$i=0,\ldots,N$$. The result becomes

(14)$\sum_{j=0}^N A^{(w)}_{i,j} c^{(w)}_j = b_i^{(w)},\quad i=0,\ldots,N,$
(15)$\sum_{j=0}^N A^{(T)}_{i,j} c^{(T)}_j = b_i^{(T)},\quad i=0,\ldots,N,$
$A^{(w)}_{i,j} = \mu(\nabla {\psi}_j,\nabla {\psi}_i),$
$b_i^{(w)} = (\beta, {\psi}_i),$
$A^{(T)}_{i,j} = \kappa(\nabla {\psi}_j,\nabla {\psi}_i),$
$b_i^{(T)} = \mu((\sum_j c^{(w)}_j\nabla{\psi}_j)\cdot (\sum_k c^{(w)}_k\nabla{\psi}_k), {\psi}_i) {\thinspace .}$

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{split}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{split}$

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

$\mu K c^{(w)} = b^{(w)},$
$\kappa K c^{(T)} = b^{(T)}{\thinspace .}$

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)}$$.

### 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

(16)$\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,$
(17)$\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,$
$A^{(w,w)}_{i,j} = \mu(\nabla {\psi}_j,{\psi}_i),$
$A^{(w,T)}_{i,j} = 0,$
$b_i^{(w)} = (\beta, {\psi}_i),$
$A^{(w,T)}_{i,j} = \mu((\nabla{\psi} w_{-})\cdot\nabla{\psi}_j), {\psi}_i),$
$A^{(T,T)}_{i,j} = \kappa(\nabla {\psi}_j,{\psi}_i),$
$b_i^{(T)} = 0 {\thinspace .}$

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

$\mu K c^{(w)} = 0 b^{(w)},$
$L c^{(w)} + \kappa K c^{(T)} =0,$

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

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

$\begin{split}\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),\end{split}$

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 (14) 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 Th 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

$\begin{split}\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 .}\end{split}$

## 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{split}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{split}$

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

(18)$\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)},$
(19)$\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)}$${\thinspace .}$

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

(20)$\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},$

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

The associated linear system is similar to (14)-(15) or (23)-(24), except that 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

(21)$\sum_{j=0}^{N_w} A^{(w)}_{i,j} c^{(w)}_j = b_i^{(w)},\quad i=0,\ldots,N_w,$
(22)$\sum_{j=0}^{N_T} A^{(T)}_{i,j} c^{(T)}_j = b_i^{(T)},\quad i=0,\ldots,N_T,$
$A^{(w)}_{i,j} = \mu(\nabla {\psi}_j^{(w)},{\psi}_i^{(w)}),$
$b_i^{(w)} = (\beta, {\psi}_i^{(w)}),$
$A^{(T)}_{i,j} = \kappa(\nabla {\psi}_j^{(T)},{\psi}_i^{(T)}),$
$b_i^{(T)} = \mu(\nabla w_{-}, {\psi}_i^{(T)}) {\thinspace .}$

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

(23)$\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,$
(24)$\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,$
$A^{(w,w)}_{i,j} = \mu(\nabla {\psi}_j^{(w)},{\psi}_i^{(w)}),$
$A^{(w,T)}_{i,j} = 0,$
$b_i^{(w)} = (\beta, {\psi}_i^{(w)}),$
$A^{(w,T)}_{i,j} = \mu (\nabla w_{-}\cdot\nabla{\psi}_j^{(w)}), {\psi}_i^{(T)}),$
$A^{(T,T)}_{i,j} = \kappa(\nabla {\psi}_j^{(T)},{\psi}_i^{(T)}),$
$b_i^{(T)} = 0 {\thinspace .}$

The corresponding block form

$\begin{split}\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),\end{split}$

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

## Computations in 1D¶

We can reduce the system (4)-(5) to one space dimension, which corresponds to flow in a channel between two flat plates. Alternatively, one may consider flow in a circular pipe, introduce cylindrical coordinates, and utilize the radial symmetry to reduce the equations to a one-dimensional problem in the radial coordinate. The former model becomes

(25)$\mu w_{xx} = -\beta,$
(26)$\kappa T_{xx} = - \mu w_x^2,$

while the model in the radial coordinate $$r$$ reads

(27)$\mu\frac{1}{r}\frac{d}{dr}\left( r\frac{dw}{dr}\right) = -\beta,$
(28)$\kappa \frac{1}{r}\frac{d}{dr}\left( r\frac{dT}{dr}\right) = - \mu \left( \frac{dw}{dr}\right)^2 {\thinspace .}$

The domain for (25)-(26) is $$\Omega = [0,H]$$, with boundary conditions $$w(0)=w(H)=0$$ and $$T(0)=T(H)=T_0$$. For (27)-(28) the domain is $$[0,R]$$ ($$R$$ being the radius of the pipe) and the boundary conditions are $$du/dr = dT/dr =0$$ for $$r=0$$, $$u(R)=0$$, and $$T(R)=T_0$$.

Calculations to be continued...

#### Previous topic

Time-dependent problems

Exercises (2)