Documentation |
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-2-by-2 tensor. By the notation $$\nabla \cdot (c\otimes \nabla u)$$, we mean the N-by-1 matrix with (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$$
The symbols a and d denote N-by-N matrices, and u denotes column vectors of lengthN.
The elements c_{ijkl}, a_{ij}, d_{ij}, and f_{i} 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 c_{ijkl} 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}}\xb7\text{\hspace{0.17em}}\left(c\otimes \nabla u\right)+qu=g+{h}^{\prime}\mu .\end{array}$$
By the notation $$n\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\left(c\otimes \nabla u\right)$$ we mean the 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)$$. 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 triangulation. Then the number of unknowns is N_{p}N = N_{u}. 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-N_{u} 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 N_{u} – M vector. The null space of H is spanned by the columns of B, and U = BV + u_{d} 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´ F – B´Ku_{d}
which is symmetric and positive definite if K is.