| Partial Differential Equation Toolbox | |
The Elliptic System
The toolbox can also handle systems of N partial differential equations over the domain 
. We have the elliptic system
- 
· (c 
u) + au = f
where c is an N-by-N-by-2-by-2 tensor. By the notation · (c u), we mean the N-by-1 matrix with (i,1)-component.
The symbols a and d denote N-by-N matrices, and u denotes column vectors 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.
The PDE Toolbox 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 procedure described above for 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,
By the notation 
we mean the N-by-1 matrix with (i,1)-component
where the outward normal vector of the boundary is 
. There are M Dirichlet conditions and the h-matrix is M-by-N, M 
0. The generalized Neumann condition contains a source h'µ 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.
The PDE Toolbox supports two implementations of Dirichlet conditions. The simplest is the "Stiff Spring" model, so named for its interpretation in solid mechanics. See the section The Elliptic Equation for the scalar case, which is equivalent to a diagonal h-matrix. For the general case, Dirichlet conditions
are approximated by adding a term
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 triangulation. 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
The constraints can be solved for m of the U-variables, the remaining called V, an Nu - m vector. The null space of H is spanned by the columns of B, and
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
which is symmetric and positive definite if K is.
| | The Elliptic Equation | The Parabolic Equation | |