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 (c\otimes \nabla u)$$ 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 (c\otimes \nabla u)$$ represents
an *N*-by-1 matrix with an (*i*,1)-component

$$\begin{array}{l}{\displaystyle \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}}\\ +{\displaystyle \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}}\\ +{\displaystyle \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 *c _{ijkl}*,

`c`

, `a`

, `d`

,
and `f`

. The case of identity, diagonal, and symmetric
matrices are handled as special cases. For the tensor 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}}\xb7\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}}\xb7\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}(\alpha ){c}_{i,j,1,1}\frac{\partial}{\partial x}+\mathrm{cos}(\alpha ){c}_{i,j,1,2}\frac{\partial}{\partial y}+\mathrm{sin}(\alpha ){c}_{i,j,2,1}\frac{\partial}{\partial x}+\mathrm{sin}(\alpha ){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}(\alpha ),\mathrm{sin}(\alpha )\right)$$.

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

$$\begin{array}{l}{\displaystyle \sum _{j=1}^{N}\left(\mathrm{cos}(\alpha ){c}_{i,j,1,1}\frac{\partial}{\partial x}+\mathrm{cos}(\alpha ){c}_{i,j,1,2}\frac{\partial}{\partial y}+\mathrm{cos}(\alpha ){c}_{i,j,1,3}\frac{\partial}{\partial z}\right)}{u}_{j}\\ +{\displaystyle \sum _{j=1}^{N}\left(\mathrm{cos}(\beta ){c}_{i,j,2,1}\frac{\partial}{\partial x}+\mathrm{cos}(\beta ){c}_{i,j,2,2}\frac{\partial}{\partial y}+\mathrm{cos}(\beta ){c}_{i,j,2,3}\frac{\partial}{\partial z}\right)}{u}_{j}\\ +{\displaystyle \sum _{j=1}^{N}\left(\mathrm{cos}(\gamma ){c}_{i,j,3,1}\frac{\partial}{\partial x}+\mathrm{cos}(\gamma ){c}_{i,j,3,2}\frac{\partial}{\partial y}+\mathrm{cos}(\gamma ){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({h}^{\prime}hu-{h}^{\prime}r)$$

to the equations **KU** = **F**, where *L* is a large number
such as 10^{4} 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 *N _{p}* nodes
in the mesh. Then the number of unknowns is

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 *N _{u}* –

*B*´ *KBV* = *B*´ *F* – *B*´*Ku _{d}*

which is symmetric and positive definite if *K* is.

Was this topic helpful?