Partial Differential Equation Toolbox™ solves equations of the form
$$m\frac{{\partial}^{2}u}{\partial {t}^{2}}+d\frac{\partial u}{\partial t}-\nabla \xb7\left(c\nabla u\right)+au=f.$$
When the m and d coefficients are 0, this reduces to
$$-\nabla \cdot \left(c\nabla u\right)+au=f,$$
which the documentation calls an elliptic equation, whether or not the equation is elliptic in the mathematical sense. The equation holds in Ω, where Ω is a bounded domain in two or three dimensions. c, a, f, and the unknown solution u are complex functions defined on Ω. c can also be a 2-by-2 matrix function on Ω. The boundary conditions specify a combination of u and its normal derivative on the boundary:
Dirichlet: hu = r on the boundary ∂Ω.
Generalized Neumann: $$\overrightarrow{n}$$ · (c∇u) + qu = g on ∂Ω.
Mixed: Only applicable to systems. A combination of Dirichlet and generalized Neumann.
$$\overrightarrow{n}$$ is the outward unit normal. g, q, h, and r are functions defined on ∂Ω.
Our nomenclature deviates slightly from the tradition for potential theory, where a Neumann condition usually refers to the case q = 0 and our Neumann would be called a mixed condition. In some contexts, the generalized Neumann boundary conditions is also referred to as the Robin boundary conditions. In variational calculus, Dirichlet conditions are also called essential boundary conditions and restrict the trial space. Neumann conditions are also called natural conditions and arise as necessary conditions for a solution. The variational form of the Partial Differential Equation Toolbox equation with Neumann conditions is given below.
The approximate solution to the elliptic PDE is found in three steps:
Describe the geometry of the domain Ω and the boundary conditions. For 2-D geometry, create geometry using the PDE app or through MATLAB^{®} files. For 3-D geometry, import the geometry in STL file format. See 2-D Geometry, Create and View 3-D Geometry, and Boundary Conditions.
Build a triangular mesh on the domain Ω. The software has mesh generating and mesh refining facilities. A mesh is described by three matrices of fixed format that contain information about the mesh points, the boundary segments, and the elements.
Discretize the PDE and the boundary conditions to obtain a linear system Ku = F. The unknown vector u contains the values of the approximate solution at the mesh points, the matrix K is assembled from the coefficients c, a, h, and q and the right-hand side F contains, essentially, averages of f around each mesh point and contributions from g. Once the matrices K and F are assembled, you have the entire MATLAB environment at your disposal to solve the linear system and further process the solution.
More elaborate applications make use of the Finite Element Method (FEM) specific information returned by the different functions of the software. Therefore we quickly summarize the theory and technique of FEM solvers to enable advanced applications to make full use of the computed quantities.
FEM can be summarized in the following sentence: Project the weak form of the differential equation onto a finite-dimensional function space. The rest of this section deals with explaining the preceding statement.
We start with the weak form of the differential equation. Without restricting the generality, we assume generalized Neumann conditions on the whole boundary, since Dirichlet conditions can be approximated by generalized Neumann conditions. In the simple case of a unit matrix h, setting g = qr and then letting q → ∞ yields the Dirichlet condition because division with a very large q cancels the normal derivative terms. The actual implementation is different, since the preceding procedure may create conditioning problems. The mixed boundary condition of the system case requires a more complicated treatment, described in Systems of PDEs.
Assume that u is a solution of the differential equation. Multiply the equation with an arbitrary test function v and integrate on Ω:
$$\underset{\Omega}{\int}\left(-\left(\nabla \text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}c\nabla u\right)v+auv\right)dx}={\displaystyle \underset{\Omega}{\int}fv\text{\hspace{0.17em}}dx}.$$
Integrate by parts (i.e., use Green's formula) to obtain
$$\underset{\Omega}{\int}\left(\left(c\nabla u\right)\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\nabla v+auv\right)\text{\hspace{0.17em}}dx}-{\displaystyle \underset{\partial \Omega}{\int}\overrightarrow{n}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\left(c\nabla u\right)v\text{\hspace{0.17em}}ds}={\displaystyle \underset{\Omega}{\int}fv\text{\hspace{0.17em}}dx}.$$
The boundary integral can be replaced by the boundary condition:
$$\underset{\Omega}{\int}\left(\left(c\nabla u\right)\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\nabla v+auv\right)\text{\hspace{0.17em}}dx}-{\displaystyle \underset{\partial \Omega}{\int}\left(-qu+g\right)v\text{\hspace{0.17em}}ds}={\displaystyle \underset{\Omega}{\int}fv\text{\hspace{0.17em}}dx}.$$
Replace the original problem with Find u such that
$$\underset{\Omega}{\int}\left(\left(c\nabla u\right)\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\nabla v+auv-fv\right)\text{\hspace{0.17em}}dx}-{\displaystyle \underset{\partial \Omega}{\int}\left(-qu+g\right)v\text{\hspace{0.17em}}ds}=0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\forall v.$$
This equation is called the variational, or weak, form of the differential equation. Obviously, any solution of the differential equation is also a solution of the variational problem. The reverse is true under some restrictions on the domain and on the coefficient functions. The solution of the variational problem is also called the weak solution of the differential equation.
The solution u and the test functions v belong to some function space V. The next step is to choose an Np-dimensional subspace $${V}_{{N}_{p}}\subset V$$. Project the weak form of the differential equation onto a finite-dimensional function space simply means requesting u and v to lie in $${V}_{{N}_{p}}$$ rather than V. The solution of the finite dimensional problem turns out to be the element of $${V}_{{N}_{p}}$$ that lies closest to the weak solution when measured in the energy norm. Convergence is guaranteed if the space $${V}_{{N}_{p}}$$ tends to V as N_{p}→∞. Since the differential operator is linear, we demand that the variational equation is satisfied for N_{p} test-functions Φ_{i }∊$${V}_{{N}_{p}}$$ that form a basis, i.e.,
$$\underset{\Omega}{\int}\left(\left(c\nabla u\right)\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\nabla {\varphi}_{i}+au{\varphi}_{i}-f{\varphi}_{i}\right)\text{\hspace{0.17em}}dx}-{\displaystyle \underset{\partial \Omega}{\int}\left(-qu+g\right){\varphi}_{i}\text{\hspace{0.17em}}ds}=0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}i=1,\mathrm{...},{N}_{p}.$$
Expand u in the same basis of $${V}_{{N}_{p}}$$ elements
$$u(x)={\displaystyle \sum _{j=1}^{{N}_{p}}{U}_{j}{\varphi}_{j}(x)},$$
and obtain the system of equations
$$\sum _{j=1}^{{N}_{p}}\left({\displaystyle \underset{\Omega}{\int}\left(\left(c\nabla {\varphi}_{j}\right)\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\nabla {\varphi}_{i}+a{\varphi}_{j}{\varphi}_{i}\right)\text{\hspace{0.17em}}dx}+{\displaystyle \underset{\partial \Omega}{\int}q{\varphi}_{j}{\varphi}_{i}\text{\hspace{0.17em}}ds}\right)}{U}_{j}={\displaystyle \underset{\Omega}{\int}f{\varphi}_{i}\text{\hspace{0.17em}}dx}+{\displaystyle \underset{\partial \Omega}{\int}g{\varphi}_{i}\text{\hspace{0.17em}}ds},\text{\hspace{0.17em}}\text{}i=1,\text{\hspace{0.17em}}\mathrm{...}\text{\hspace{0.17em}},{N}_{p.$$
Use the following notations:
$${K}_{i,j}={\displaystyle \underset{\Omega}{\int}\left(c\nabla {\varphi}_{j}\right)\text{\hspace{0.17em}}\cdot \text{\hspace{0.17em}}\nabla {\varphi}_{i}\text{\hspace{0.17em}}dx}\text{\hspace{1em}}\text{(stiffnessmatrix)}$$
$${M}_{i,j}={\displaystyle \underset{\Omega}{\int}a{\varphi}_{j}{\varphi}_{i}\text{\hspace{0.17em}}dx}\text{\hspace{1em}}\text{(massmatrix)}$$
$${Q}_{i,j}={\displaystyle \underset{\partial \Omega}{\int}q{\varphi}_{j}{\varphi}_{i}\text{\hspace{0.17em}}ds}$$
$${F}_{i}={\displaystyle \underset{\Omega}{\int}f{\varphi}_{i}\text{\hspace{0.17em}}dx}$$
$${G}_{i}={\displaystyle \underset{\partial \Omega}{\int}g{\varphi}_{i}\text{\hspace{0.17em}}ds}$$
and rewrite the system in the form
(K + M + Q)U = F + G.
K, M, and Q are N_{p}-by-N_{p} matrices,
and F and G are N_{p}-vectors. K, M,
and F are produced by assema
,
while Q, G are produced by assemb
.
When it is not necessary to distinguish K, M,
and Q or F and G,
we collapse the notations to KU = F,
which form the output of assempde
.
When the problem is self-adjoint and elliptic in the usual mathematical sense, the matrix K + M + Q becomes symmetric and positive definite. Many common problems have these characteristics, most notably those that can also be formulated as minimization problems. For the case of a scalar equation, K, M, and Q are obviously symmetric. If c(x) ≥ δ > 0, a(x) ≥ 0 and q(x) ≥ 0 with q(x) > 0 on some part of ∂Ω, then, if U ≠ 0.
$${U}^{T}\left(K+M+Q\right)U={\displaystyle \underset{\Omega}{\int}\left(c{\left|u\right|}^{2}+a{u}^{2}\right)\text{\hspace{0.17em}}dx}+{\displaystyle \underset{\partial \Omega}{\int}q{u}^{2}\text{\hspace{0.17em}}ds}>0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{if}U\ne 0.$$
U^{T}(K + M + Q)U is the energy norm. There are many choices of the test-function spaces. The software uses continuous functions that are linear on each element of a 2-D mesh, and are linear or quadratic on elements of a 3-D mesh. Piecewise linearity guarantees that the integrals defining the stiffness matrix K exist. Projection onto $${V}_{{N}_{p}}$$ is nothing more than linear interpolation, and the evaluation of the solution inside an element is done just in terms of the nodal values. If the mesh is uniformly refined, $${V}_{{N}_{p}}$$ approximates the set of smooth functions on Ω.
A suitable basis for $${V}_{{N}_{p}}$$ in 2-D is the set of "tent" or "hat" functions ϕ_{i}. These are linear on each element and take the value 0 at all nodes x_{j} except for x_{i}. For the definition of basis functions for 3-D geometry, see Finite Element Basis for 3-D. Requesting ϕ_{i}(x_{i}) = 1 yields the very pleasant property
$$u\left({x}_{i}\right)={\displaystyle \sum _{j=1}^{{N}_{p}}{U}_{j}{\varphi}_{j}\left({x}_{i}\right)}={U}_{i}.$$
That is, by solving the FEM system we obtain the nodal values of the approximate solution. The basis function ϕ_{i} vanishes on all the elements that do not contain the node x_{i}. The immediate consequence is that the integrals appearing in K_{i,j}, M_{i,j}, Q_{i,j}, F_{i} and G_{i} only need to be computed on the elements that contain the node x_{i}. Secondly, it means that K_{i,j} andM_{i,j} are zero unless x_{i} and x_{j} are vertices of the same element and thus K and M are very sparse matrices. Their sparse structure depends on the ordering of the indices of the mesh points.
The integrals in the FEM matrices are computed by adding the
contributions from each element to the corresponding entries (i.e.,
only if the corresponding mesh point is a vertex of the element).
This process is commonly called assembling, hence
the name of the function assempde
.
The assembling routines scan the elements of the mesh. For each element they compute the so-called local matrices and add their components to the correct positions in the sparse matrices or vectors.
The discussion now specializes to triangular meshes in 2-D. The local 3-by-3 matrices contain the integrals evaluated only on the current triangle. The coefficients are assumed constant on the triangle and they are evaluated only in the triangle barycenter. The integrals are computed using the midpoint rule. This approximation is optimal since it has the same order of accuracy as the piecewise linear interpolation.
Consider a triangle given by the nodes P_{1}, P_{2}, and P_{3} as in the following figure.
The Local Triangle P1P2P3
Note: The local 3-by-3 matrices contain the integrals evaluated only on the current triangle. The coefficients are assumed constant on the triangle and they are evaluated only in the triangle barycenter. |
The simplest computations are for the local mass matrix m:
$${m}_{i,j}={\displaystyle \underset{\Delta {P}_{1}{P}_{2}{P}_{3}}{\int}a\left({P}_{c}\right){\varphi}_{i}\left(x\right){\varphi}_{j}\left(x\right)\text{\hspace{0.17em}}dx}=a\left({P}_{c}\right)\frac{\text{area}\left(\Delta {P}_{1}{P}_{2}{P}_{3}\right)}{12}\left(1+{\delta}_{i,j}\right),$$
where P_{c} is the center of mass of Δ P_{1}P_{2}P_{3}, i.e.,
$${P}_{c}=\frac{{P}_{1}+{P}_{2}+{P}_{3}}{3}.$$
The contribution to the right side F is just
$${f}_{i}=f\left({P}_{c}\right)\frac{\text{area}\left(\Delta {P}_{1}{P}_{2}{P}_{3}\right)}{3}.$$
For the local stiffness matrix we have to evaluate the gradients of the basis functions that do not vanish on P_{1}P_{2}P_{3}. Since the basis functions are linear on the triangle P_{1}P_{2}P_{3}, the gradients are constants. Denote the basis functions ϕ_{1}, ϕ_{2}, and ϕ_{3} such that ϕ(P_{i}) = 1. If P_{2} – P_{3} = [x_{1},y_{1}]^{T} then we have that
$$\nabla {\varphi}_{1}=\frac{1}{2\text{\hspace{0.17em}}\text{area}\left(\Delta {P}_{1}{P}_{2}{P}_{3}\right)}\left[\begin{array}{c}{y}_{1}\\ -{x}_{1}\end{array}\right]$$
and after integration (taking c as a constant matrix on the triangle)
$${k}_{i,j}=\frac{1}{4\text{\hspace{0.17em}}\text{area}\left(\Delta {P}_{1}{P}_{2}{P}_{3}\right)}\left[{y}_{j},-{x}_{j}\right]c\left({P}_{c}\right)\left[\begin{array}{c}{y}_{1}\\ -{x}_{1}\end{array}\right].$$
If two vertices of the triangle lie on the boundary ∂Ω, they contribute to the line integrals associated to the boundary conditions. If the two boundary points are P_{1} and P_{2}, then we have
$${Q}_{i,j}=q\left({P}_{b}\right)\frac{\Vert {P}_{1}-{P}_{2}\Vert}{6}\left(1+{\delta}_{i,j}\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}i,j=1,2$$
and
$${G}_{i}=g\left({P}_{b}\right)\frac{\Vert {P}_{1}-{P}_{2}\Vert}{2},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}i=1,2$$
where P_{b} is the midpoint of P_{1}P_{2}.
For each triangle the vertices P_{m} of the local triangle correspond to the indices i_{m} of the mesh points. The contributions of the individual triangle are added to the matrices such that, e.g.,
$${K}_{{i}_{m},{i}_{n}}t\leftarrow {K}_{{i}_{m},{i}_{n}}+{k}_{m,n},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}m,n=1,2,3.$$
This is done by the function assempde
.
The gradients and the areas of the triangles are computed by the function pdetrg
.
The Dirichlet boundary conditions are treated in a slightly
different manner. They are eliminated from the linear system by a
procedure that yields a symmetric, reduced system. The function assempde
can
return matrices K, F, B,
and ud such that the solution is u = Bv + ud where Kv = F. u is
an N_{p}-vector, and if the
rank of the Dirichlet conditions is rD, then v has N_{p} – rD components.