## Documentation Center |

The basic elliptic equation handled by the software is

in Ω, where Ω is a bounded domain in the plane. *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:*· (*c*∇*u*) +*qu*=*g*on ∂Ω.*Mixed:*Only applicable to*systems*. A combination of Dirichlet and generalized Neumann.

*
* 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. This can be done either interactively using the PDE app or through MATLAB

^{®}files (see`pdegeom`and`pdebound`).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 triangles.

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 Ω:

Integrate by parts (i.e., use Green's formula) to obtain

The boundary integral can be replaced by the boundary condition:

Replace the original problem with *Find u such that*

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
. *Project
the weak form of the differential equation onto a finite-dimensional
function space* simply means requesting *u* and *v* to
lie in
rather than *V*.
The solution of the finite dimensional problem turns out to be the
element of
that lies closest to the weak
solution when measured in the energy norm. Convergence is guaranteed
if the space
tends to *V* as *N _{p}*→∞.
Since the differential operator is linear, we demand that the variational
equation is satisfied for

Expand *u* in the same basis of
elements

and obtain the system of equations

Use the following notations:

and rewrite the system in the form

(*K* + *M* + *Q*)*U* = *F* + *G*.

*K*, *M*, and *Q* are *N _{p}*-by-

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}*(

A suitable basis for
is
the set of "tent" or "hat" functions *ϕ*_{i}.
These are linear on each triangle and take the value 0 at all nodes *x _{j}* except
for

That is, by solving the FEM system we obtain the nodal values
of the approximate solution. The basis function *ϕ*_{i} vanishes
on all the triangles that do not contain the node *x _{i}*.
The immediate consequence is that the integrals appearing in

The integrals in the FEM matrices are computed by adding the
contributions from each triangle to the corresponding entries (i.e.,
only if the corresponding mesh point is a vertex of the triangle).
This process is commonly called *assembling*, hence
the name of the function `assempde`.

The assembling routines scan the triangles of the mesh. For each triangle they compute the so-called local matrices and add their components to the correct positions in the sparse matrices or vectors. (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 simplest computations are for the local mass matrix *m*:

where *P _{c}* is the center
of mass of Δ

The contribution to the right side *F* is just

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

and after integration (taking *c* as a constant
matrix on the triangle)

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

and

where *P _{b}* is the midpoint
of

For each triangle the vertices *P _{m}* of
the local triangle correspond to the indices

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

Was this topic helpful?