Partial Differential Equation Toolbox™

The Elliptic System

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

the parabolic system

the hyperbolic system

and the eigenvalue system

where is an N-by-N-by-2-by-2 tensor. By the notation , we mean the N-by-1 matrix with (i,1)-component.

The symbols and denote N-by-N matrices, and denotes column vectors of length N.

The elements c_ijkl, a_ij, d_ij, and f_i of , , , and 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,

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 -matrix is M-by-N, M ≥ 0. The generalized Neumann condition contains a source , 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 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 10⁴ times a representative size of the elements of K.

When this number is increased, 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 , and is free of the ill-conditioning, but is more involved computationally. Assume that there are nodes in the triangulation. Then the number of unknowns is N_pN = 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, . 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.

The Elliptic Equation

The Parabolic Equation