# Documentation

## Systems of PDEs

Partial Differential Equation Toolbox™ software can also handle systems of N partial differential equations over the domain Ω. We have the elliptic system

$-\nabla \cdot \left(c\otimes \nabla u\right)+au=f.$

the parabolic system

$d\frac{\partial u}{\partial t}-\nabla \cdot \left(c\otimes \nabla u\right)+au=f,$

the hyperbolic system

$d\frac{{\partial }^{2}u}{\partial {t}^{2}}-\nabla \cdot \left(c\otimes \nabla u\right)+au=f,$

and the eigenvalue system

$-\nabla \cdot \left(c\otimes \nabla u\right)+au=\lambda du,$

where c is an N-by-N-by-D-by-D tensor, and D is the geometry dimensions, 2 or 3.

For 2-D systems, the notation $\nabla \cdot \left(c\otimes \nabla u\right)$ represents an N-by-1 matrix with an (i,1)-component

$\sum _{j=1}^{N}\left(\frac{\partial }{\partial x}{c}_{i,j,1,1}\frac{\partial }{\partial x}+\frac{\partial }{\partial x}{c}_{i,j,1,2}\frac{\partial }{\partial y}+\frac{\partial }{\partial y}{c}_{i,j,2,1}\frac{\partial }{\partial x}+\frac{\partial }{\partial y}{c}_{i,j,2,2}\frac{\partial }{\partial y}\right){u}_{j}.$

For 3-D systems, the notation $\nabla \cdot \left(c\otimes \nabla u\right)$ represents an N-by-1 matrix with an (i,1)-component

$\begin{array}{l}\sum _{j=1}^{N}\left(\frac{\partial }{\partial x}{c}_{i,j,1,1}\frac{\partial }{\partial x}+\frac{\partial }{\partial x}{c}_{i,j,1,2}\frac{\partial }{\partial y}+\frac{\partial }{\partial x}{c}_{i,j,1,3}\frac{\partial }{\partial z}\right){u}_{j}\\ +\sum _{j=1}^{N}\left(\frac{\partial }{\partial y}{c}_{i,j,2,1}\frac{\partial }{\partial x}+\frac{\partial }{\partial y}{c}_{i,j,2,2}\frac{\partial }{\partial y}+\frac{\partial }{\partial y}{c}_{i,j,2,3}\frac{\partial }{\partial z}\right){u}_{j}\\ +\sum _{j=1}^{N}\left(\frac{\partial }{\partial z}{c}_{i,j,3,1}\frac{\partial }{\partial x}+\frac{\partial }{\partial z}{c}_{i,j,3,2}\frac{\partial }{\partial y}+\frac{\partial }{\partial z}{c}_{i,j,3,3}\frac{\partial }{\partial z}\right){u}_{j}.\end{array}$

The symbols a and d denote N-by-N matrices, and f denotes a column vector of length N.

The elements cijkl, aij, dij, and fi of c, a, d, and f are stored row-wise in the MATLAB® matrices `c`, `a`, `d`, and `f`. The case of identity, diagonal, and symmetric matrices are handled as special cases. For the tensor cijkl this applies both to the indices i and j, and to the indices k and l.

Partial Differential Equation Toolbox software does not check the ellipticity of the problem, and it is quite possible to define a system that is not elliptic in the mathematical sense. The preceding procedure that describes the scalar case is applied to each component of the system, yielding a symmetric positive definite system of equations whenever the differential system possesses these characteristics.

The boundary conditions now in general are mixed, i.e., for each point on the boundary a combination of Dirichlet and generalized Neumann conditions,

$\begin{array}{l}hu=r\\ n\text{\hspace{0.17em}}·\text{\hspace{0.17em}}\left(c\otimes \nabla u\right)+qu=g+{h}^{\prime }\mu .\end{array}$

For 2-D systems, the notation $n\text{\hspace{0.17em}}·\text{\hspace{0.17em}}\left(c\otimes \nabla u\right)$ represents an N-by-1 matrix with (i,1)-component

$\sum _{j=1}^{N}\left(\mathrm{cos}\left(\alpha \right){c}_{i,j,1,1}\frac{\partial }{\partial x}+\mathrm{cos}\left(\alpha \right){c}_{i,j,1,2}\frac{\partial }{\partial y}+\mathrm{sin}\left(\alpha \right){c}_{i,j,2,1}\frac{\partial }{\partial x}+\mathrm{sin}\left(\alpha \right){c}_{i,j,2,2}\frac{\partial }{\partial y}\right){u}_{j}$

where the outward normal vector of the boundary is $n=\left(\mathrm{cos}\left(\alpha \right),\mathrm{sin}\left(\alpha \right)\right)$.

For 3-D systems, the notation $n\text{\hspace{0.17em}}·\text{\hspace{0.17em}}\left(c\otimes \nabla u\right)$ represents an N-by-1 matrix with (i,1)-component

$\begin{array}{l}\sum _{j=1}^{N}\left(\mathrm{cos}\left(\alpha \right){c}_{i,j,1,1}\frac{\partial }{\partial x}+\mathrm{cos}\left(\alpha \right){c}_{i,j,1,2}\frac{\partial }{\partial y}+\mathrm{cos}\left(\alpha \right){c}_{i,j,1,3}\frac{\partial }{\partial z}\right){u}_{j}\\ +\sum _{j=1}^{N}\left(\mathrm{cos}\left(\beta \right){c}_{i,j,2,1}\frac{\partial }{\partial x}+\mathrm{cos}\left(\beta \right){c}_{i,j,2,2}\frac{\partial }{\partial y}+\mathrm{cos}\left(\beta \right){c}_{i,j,2,3}\frac{\partial }{\partial z}\right){u}_{j}\\ +\sum _{j=1}^{N}\left(\mathrm{cos}\left(\gamma \right){c}_{i,j,3,1}\frac{\partial }{\partial x}+\mathrm{cos}\left(\gamma \right){c}_{i,j,3,2}\frac{\partial }{\partial y}+\mathrm{cos}\left(\gamma \right){c}_{i,j,3,3}\frac{\partial }{\partial z}\right){u}_{j},\end{array}$

where the outward normal to the boundary is

$n=\left(\mathrm{cos}\left(\alpha \right),\mathrm{cos}\left(\beta \right),\mathrm{cos}\left(\gamma \right)\right).$

There are M Dirichlet conditions and the h-matrix is M-by-N, M ≥ 0. The generalized Neumann condition contains a source ${h}^{\prime }\mu$, where the Lagrange multipliers μ are computed such that the Dirichlet conditions become satisfied. In a structural mechanics problem, this term is exactly the reaction force necessary to satisfy the kinematic constraints described by the Dirichlet conditions.

The rest of this section details the treatment of the Dirichlet conditions and may be skipped on a first reading.

Partial Differential Equation Toolbox software supports two implementations of Dirichlet conditions. The simplest is the "Stiff Spring" model, so named for its interpretation in solid mechanics. See Elliptic Equations for the scalar case, which is equivalent to a diagonal h-matrix. For the general case, Dirichlet conditions

hu = r

are approximated by adding a term

$L\left({h}^{\prime }hu-{h}^{\prime }r\right)$

to the equations KU = F, where L is a large number such as 104 times a representative size of the elements of K.

When this number is increased, hu = r will be more accurately satisfied, but the potential ill-conditioning of the modified equations will become more serious.

The second method is also applicable to general mixed conditions with nondiagonal h, and is free of the ill-conditioning, but is more involved computationally. Assume that there are Np nodes in the mesh. Then the number of unknowns is NpN = Nu. When Dirichlet boundary conditions fix some of the unknowns, the linear system can be correspondingly reduced. This is easily done by removing rows and columns when u values are given, but here we must treat the case when some linear combinations of the components of u are given, hu = r. These are collected into HU = R where H is an M-by-Nu matrix and R is an M-vector.

With the reaction force term the system becomes

KU +H´ µ = F

HU = R.

The constraints can be solved for M of the U-variables, the remaining called V, an NuM vector. The null space of H is spanned by the columns of B, and U = BV + ud makes U satisfy the Dirichlet conditions. A permutation to block-diagonal form exploits the sparsity of H to speed up the following computation to find B in a numerically stable way. µ can be eliminated by premultiplying by B´ since, by the construction, HB = 0 or B´H´ = 0. The reduced system becomes

B´ KBV = B´ FB´Kud

which is symmetric and positive definite if K is.