(Not recommended) Assemble finite element matrices and solve elliptic PDE

`assempde`

is not recommended. Use `solvepde`

instead.

solves the PDE`u`

= assempde(`model`

,`c`

,`a`

,`f`

)

$$-\nabla \cdot \left(c\nabla u\right)+au=f$$

with geometry, boundary conditions, and finite element mesh in
`model`

, and coefficients `c`

,
`a`

, and `f`

. If the PDE is a system of
equations (`model.PDESystemSize`

> 1), then
`assempde`

solves the system of equations

$$-\nabla \cdot \left(c\otimes \nabla u\right)+au=f$$

`[`

,
for any of the previous input syntaxes, assembles finite element matrices using
the `Kc`

,`Fc`

,`B`

,`ud`

]
= assempde(___)*reduced linear system* form, which eliminates any
Dirichlet boundary conditions from the system of linear equations. You can
calculate the solution `u`

at node points by the command
`u = B*(Kc\Fc) + ud`

. See Reduced Linear System.

`[`

assembles finite element matrices that represent any Dirichlet boundary
conditions using a `Ks`

,`Fs`

]
= assempde(___)*stiff-spring* approximation. You can
calculate the solution `u`

at node points by the command
`u = Ks\Fs`

. See Stiff-Spring Approximation.

Solve an elliptic PDE on an L-shaped region.

Create a scalar PDE model. Incorporate the geometry of an L-shaped region.

model = createpde; geometryFromEdges(model,@lshapeg);

Apply zero Dirichlet boundary conditions to all edges.

applyBoundaryCondition(model,'Edge',1:model.Geometry.NumEdges,'u',0);

Generate a finite element mesh.

generateMesh(model,'GeometricOrder','linear');

Solve the PDE $$-\nabla \cdot (c\nabla u)+au=f$$ with parameters `c = 1`

, `a = 0`

, and `f = 5`

.

c = 1; a = 0; f = 5; u = assempde(model,c,a,f);

Plot the solution.

`pdeplot(model,'XYData',u)`

Solve a 3-D elliptic PDE using a PDE model.

Create a PDE model container, import a 3-D geometry description, and view the geometry.

model = createpde; importGeometry(model,'Block.stl'); pdegplot(model,'FaceLabels','on','FaceAlpha',0.5)

Set zero Dirichlet conditions on faces 1 through 4 (the edges). Set Neumann conditions with `g = -1`

on face 6 and `g = 1`

on face 5.

applyBoundaryCondition(model,'Face',1:4,'u',0); applyBoundaryCondition(model,'Face',6,'g',-1); applyBoundaryCondition(model,'Face',5,'g',1);

Set coefficients `c = 1`

, `a = 0`

, and `f = 0.1`

.

c = 1; a = 0; f = 0.1;

Create a mesh and solve the problem.

generateMesh(model); u = assempde(model,c,a,f);

Plot the solution on the surface.

`pdeplot3D(model,'ColorMapData',u)`

Solve a 2-D PDE using the older syntax for mesh.

Create a circle geometry.

g = @circleg;

Set zero Dirichlet boundary conditions.

b = @circleb1;

Create a mesh for the geometry.

[p,e,t] = initmesh(g);

Solve the PDE $$-\nabla \cdot (c\nabla u)+au=f$$ with parameters `c = 1`

, `a = 0`

, and `f = sin(x)`

.

```
c = 1;
a = 0;
f = 'sin(x)';
u = assempde(b,p,e,t,c,a,f);
```

Plot the solution.

`pdeplot(p,e,t,'XYData',u)`

Obtain the finite-element matrices that represent the problem using a reduced linear algebra representation of Dirichlet boundary conditions.

Create a scalar PDE model. Import a simple 3-D geometry.

```
model = createpde;
importGeometry(model,'Block.stl');
```

Set zero Dirichlet boundary conditions on all the geometry faces.

applyBoundaryCondition(model,'dirichlet','Face',1:model.Geometry.NumFaces,'u',0);

Generate a mesh for the geometry.

generateMesh(model);

Obtain finite element matrices `K`

, `F`

, `B`

, and `ud`

that represent the equation $$-\nabla \cdot (c\nabla u)+au=f$$ with parameters $$c=1$$, $$a=0$$, and $$f=\mathrm{log}(1+x+\frac{y}{1+z})$$.

```
c = 1;
a = 0;
f = 'log(1+x+y./(1+z))';
[K,F,B,ud] = assempde(model,c,a,f);
```

You can obtain the solution `u`

of the PDE at mesh nodes by executing the command

u = B*(K\F) + ud;

Generally, this solution is slightly more accurate than the stiff-spring solution, as calculated in the next example.

Obtain the stiff-spring approximation of finite element matrices.

Create a scalar PDE model. Import a simple 3-D geometry.

```
model = createpde;
importGeometry(model,'Block.stl');
```

Set zero Dirichlet boundary conditions on all the geometry faces.

applyBoundaryCondition(model,'Face',1:model.Geometry.NumFaces,'u',0);

Generate a mesh for the geometry.

generateMesh(model);

Obtain finite element matrices `Ks`

and `Fs`

that represent the equation $$-\nabla \cdot (c\nabla u)+au=f$$ with parameters $$c=1$$, $$a=0$$, and $$f=\mathrm{log}(1+x+\frac{y}{1+z})$$.

```
c = 1;
a = 0;
f = 'log(1+x+y./(1+z))';
[Ks,Fs] = assempde(model,c,a,f);
```

You can obtain the solution `u`

of the PDE at mesh nodes by executing the command

u = Ks\Fs;

Generally, this solution is slightly less accurate than the reduced linear algebra solution, as calculated in the previous example.

Obtain the full collection of finite element matrices for an elliptic problem.

Import geometry and set up an elliptic problem with Dirichlet boundary conditions. The `Torus.stl`

geometry has only one face, so you need set only one boundary condition.

model = createpde(); importGeometry(model,'Torus.stl'); applyBoundaryCondition(model,'face',1,'u',0); c = 1; a = 0; f = 1; generateMesh(model);

Create the finite element matrices that represent this problem.

[K,M,F,Q,G,H,R] = assempde(model,c,a,f);

Most of the resulting matrices are quite sparse. `G`

, `M`

, `Q`

, and `R`

are all zero sparse matrices.

howsparse = @(x)nnz(x)/numel(x); disp(['Maximum fraction of nonzero entries in K or H is ',... num2str(max(howsparse(K),howsparse(H)))])

Maximum fraction of nonzero entries in K or H is 0.002006

To find the solution to the PDE, call `assempde`

again.

u = assempde(K,M,F,Q,G,H,R);

`model`

— PDE model`PDEModel`

objectPDE model, specified as a `PDEModel`

object.

**Example: **`model = createpde`

`c`

— PDE coefficientscalar | matrix | character vector | character array | string scalar | string vector | coefficient function

PDE coefficient, specified as a scalar, matrix, character vector, character array, string
scalar, string vector, or coefficient function. `c`

represents the *c* coefficient in the scalar PDE

$$-\nabla \cdot \left(c\nabla u\right)+au=f$$

or in the system of PDEs

$$-\nabla \cdot \left(c\otimes \nabla u\right)+au=f$$

**Example: **`'cosh(x+y.^2)'`

**Data Types: **`double`

| `char`

| `string`

| `function_handle`

**Complex Number Support: **Yes

`a`

— PDE coefficientscalar | matrix | character vector | character array | string scalar | string vector | coefficient function

PDE coefficient, specified as a scalar, matrix, character vector, character array, string
scalar, string vector, or coefficient function. `a`

represents the
*a* coefficient in the scalar PDE

$$-\nabla \cdot \left(c\nabla u\right)+au=f$$

or in the system of PDEs

$$-\nabla \cdot \left(c\otimes \nabla u\right)+au=f$$

**Example: **`2*eye(3)`

**Data Types: **`double`

| `char`

| `string`

| `function_handle`

**Complex Number Support: **Yes

`f`

— PDE coefficientscalar | matrix | character vector | character array | string scalar | string vector | coefficient function

PDE coefficient, specified as a scalar, matrix, character vector, character array, string
scalar, string vector, or coefficient function. `f`

represents the *f* coefficient in the scalar PDE

$$-\nabla \cdot \left(c\nabla u\right)+au=f$$

or in the system of PDEs

$$-\nabla \cdot \left(c\otimes \nabla u\right)+au=f$$

**Example: **`char('sin(x)';'cos(y)';'tan(z)')`

**Data Types: **`double`

| `char`

| `string`

| `function_handle`

**Complex Number Support: **Yes

`b`

— Boundary conditionsboundary matrix | boundary file

Boundary conditions, specified as a boundary matrix or boundary file. Pass a boundary file as a function handle or as a file name. A boundary matrix is generally an export from the PDE Modeler app.

**Example: **`b = 'circleb1'`

, `b = "circleb1"`

, or ```
b =
@circleb1
```

**Data Types: **`double`

| `char`

| `string`

| `function_handle`

`p`

— Mesh pointsmatrix

Mesh points, specified as a 2-by-`Np`

matrix
of points, where `Np`

is the number of points in
the mesh. For a description of the (`p`

,`e`

,`t`

)
matrices, see Mesh Data.

Typically, you use the `p`

, `e`

, and `t`

data exported from the **PDE Modeler** app, or generated by `initmesh`

or `refinemesh`

.

**Example: **`[p,e,t] = initmesh(gd)`

**Data Types: **`double`

`e`

— Mesh edgesmatrix

Mesh edges, specified as a `7`

-by-`Ne`

matrix
of edges, where `Ne`

is the number of edges in the
mesh. For a description of the (`p`

,`e`

,`t`

)
matrices, see Mesh Data.

Typically, you use the `p`

, `e`

, and `t`

data exported from the **PDE Modeler** app, or generated by `initmesh`

or `refinemesh`

.

**Example: **`[p,e,t] = initmesh(gd)`

**Data Types: **`double`

`t`

— Mesh trianglesmatrix

Mesh triangles, specified as a `4`

-by-`Nt`

matrix
of triangles, where `Nt`

is the number of triangles
in the mesh. For a description of the (`p`

,`e`

,`t`

)
matrices, see Mesh Data.

Typically, you use the `p`

, `e`

, and `t`

data exported from the **PDE Modeler** app, or generated by `initmesh`

or `refinemesh`

.

**Example: **`[p,e,t] = initmesh(gd)`

**Data Types: **`double`

`K`

— Stiffness matrixsparse matrix | full matrix

Stiffness matrix, specified as a sparse matrix or full matrix. Generally,
you obtain `K`

from a previous call to `assema`

or
`assempde`

. For the meaning of stiffness matrix, see
Elliptic Equations.

**Example: **```
[K,M,F,Q,G,H,R] =
assempde(model,c,a,f)
```

**Data Types: **`double`

**Complex Number Support: **Yes

`M`

— Mass matrixsparse matrix | full matrix

Mass matrix, specified as a sparse matrix or full matrix. Generally, you
obtain `M`

from a previous call to `assema`

or
`assempde`

. For the meaning of mass matrix, see Elliptic Equations.

**Example: **```
[K,M,F,Q,G,H,R] =
assempde(model,c,a,f)
```

**Data Types: **`double`

**Complex Number Support: **Yes

`F`

— Finite element `f`

representationvector

Finite element `f`

representation, specified as a vector.
Generally, you obtain `F`

from a previous call to `assema`

or
`assempde`

. For the meaning of this representation,
see Elliptic Equations.

**Example: **```
[K,M,F,Q,G,H,R] =
assempde(model,c,a,f)
```

**Data Types: **`double`

**Complex Number Support: **Yes

`Q`

— Neumann boundary condition matrixsparse matrix | full matrix

Neumann boundary condition matrix, specified as a sparse matrix or full
matrix. Generally, you obtain `Q`

from a previous call to
`assemb`

or
`assempde`

. For the meaning of this matrix, see Elliptic Equations.

**Example: **```
[K,M,F,Q,G,H,R] =
assempde(model,c,a,f)
```

**Data Types: **`double`

**Complex Number Support: **Yes

`G`

— Neumann boundary condition vectorsparse vector | full vector

Neumann boundary condition vector, specified as a sparse vector or full
vector. Generally, you obtain `G`

from a previous call to
`assemb`

or
`assempde`

. For the meaning of this vector, see Elliptic Equations.

**Example: **```
[K,M,F,Q,G,H,R] =
assempde(model,c,a,f)
```

**Data Types: **`double`

**Complex Number Support: **Yes

`H`

— Dirichlet boundary condition matrixsparse matrix | full matrix

Dirichlet boundary condition matrix, specified as a sparse matrix or full
matrix. Generally, you obtain `H`

from a previous call to
`assemb`

or
`assempde`

. For the meaning of this matrix, see Algorithms.

**Example: **```
[K,M,F,Q,G,H,R] =
assempde(model,c,a,f)
```

**Data Types: **`double`

**Complex Number Support: **Yes

`R`

— Dirichlet boundary condition vectorsparse vector | full vector

Dirichlet boundary condition vector, specified as a sparse vector or full
vector. Generally, you obtain `R`

from a previous call to
`assemb`

or
`assempde`

. For the meaning of this vector, see Algorithms.

**Example: **```
[K,M,F,Q,G,H,R] =
assempde(model,c,a,f)
```

**Data Types: **`double`

**Complex Number Support: **Yes

`sdl`

— Subdomain labelsvector of positive integers

Subdomain labels, specified as a vector of positive integers. For 2-D geometry only. View the subdomain labels in your geometry using the command

pdegplot(g,'SubdomainLabels','on')

**Example: **`sdl = [1,3:5];`

**Data Types: **`double`

`u`

— PDE solutionvector

PDE solution, returned as a vector.

If the PDE is scalar, meaning only one equation, then

`u`

is a column vector representing the solution*u*at each node in the mesh.`u(i)`

is the solution at the`i`

th column of`model.Mesh.Nodes`

or the`i`

th column of`p`

.If the PDE is a system of

*N*> 1 equations, then`u`

is a column vector with*N**`Np`

elements, where`Np`

is the number of nodes in the mesh. The first`Np`

elements of`u`

represent the solution of equation 1, then next`Np`

elements represent the solution of equation 2, etc.

To obtain the solution at an arbitrary point in the geometry,
use `pdeInterpolant`

.

To plot the solution, use `pdeplot`

for
2-D geometry, or see Plot 3-D Solutions and Their Gradients.

`Kc`

— Stiffness matrixsparse matrix

Stiffness matrix, returned as a sparse matrix. See Elliptic Equations.

`u1 = Kc\Fc`

returns the solution on the non-Dirichlet
points. To obtain the solution `u`

at the nodes of the
mesh,

`u = B*(Kc\Fc) + ud`

Generally, `Kc`

, `Fc`

,
`B`

, and `ud`

make a slower but more
accurate solution than `Ks`

and
`Fs`

.

`Fc`

— Load vectorvector

Load vector, returned as a vector. See Elliptic Equations.

`u = B*(Kc\Fc) + ud`

Generally, `Kc`

, `Fc`

,
`B`

, and `ud`

make a slower but more
accurate solution than `Ks`

and
`Fs`

.

`B`

— Dirichlet nullspacesparse matrix

Dirichlet nullspace, returned as a sparse matrix. See Algorithms.

`u = B*(Kc\Fc) + ud`

Generally, `Kc`

, `Fc`

,
`B`

, and `ud`

make a slower but more
accurate solution than `Ks`

and
`Fs`

.

`ud`

— Dirichlet vectorvector

Dirichlet vector, returned as a vector. See Algorithms.

`u = B*(Kc\Fc) + ud`

Generally, `Kc`

, `Fc`

,
`B`

, and `ud`

make a slower but more
accurate solution than `Ks`

and
`Fs`

.

`Ks`

— Stiffness matrix corresponding to the stiff-spring approximation for Dirichlet boundary conditionsparse matrix

Finite element matrix for stiff-spring approximation, returned as a sparse matrix. See Algorithms.

To obtain the solution `u`

at the nodes of the
mesh,

`u = Ks\Fs`

.

Generally, `Ks`

and `Fs`

make a quicker
but less accurate solution than `Kc`

,
`Fc`

, `B`

, and
`ud`

.

`Fs`

— Load vector corresponding to the stiff-spring approximation for Dirichlet boundary conditionvector

Load vector corresponding to the stiff-spring approximation for Dirichlet boundary condition, returned as a vector. See Algorithms.

To obtain the solution `u`

at the nodes of the
mesh,

`u = Ks\Fs`

.

Generally, `Ks`

and `Fs`

make a quicker
but less accurate solution than `Kc`

,
`Fc`

, `B`

, and
`ud`

.

`K`

— Stiffness matrixsparse matrix

Stiffness matrix, returned as a sparse matrix. See Elliptic Equations.

`K`

represents the stiffness matrix alone, unlike
`Kc`

or `Ks`

, which are stiffness
matrices combined with other terms to enable immediate solution of a
PDE.

Typically, you use `K`

in a subsequent call to a solver
such as `assempde`

or
`hyperbolic`

.

`M`

— Mass matrixsparse matrix

Mass matrix. returned as a sparse matrix. See Elliptic Equations.

Typically, you use `M`

in a subsequent call
to a solver such as `assempde`

or `hyperbolic`

.

`F`

— Load vectorvector

Load vector, returned as a vector. See Elliptic Equations.

`F`

represents the load vector alone, unlike `Fc`

or `Fs`

, which are load vectors
combined with other terms to enable immediate solution of a PDE.

Typically, you use `F`

in a subsequent call to a solver
such as `assempde`

or
`hyperbolic`

.

`Q`

— Neumann boundary condition matrixsparse matrix

Neumann boundary condition matrix, returned as a sparse matrix. See Elliptic Equations.

Typically, you use `Q`

in a subsequent call
to a solver such as `assempde`

or `hyperbolic`

.

`G`

— Neumann boundary condition vectorsparse vector

Neumann boundary condition vector, returned as a sparse vector. See Elliptic Equations.

Typically, you use `G`

in a subsequent call
to a solver such as `assempde`

or `hyperbolic`

.

`H`

— Dirichlet matrixsparse matrix

Dirichlet matrix, returned as a sparse matrix. See Algorithms.

Typically, you use `H`

in a subsequent call to a solver
such as `assempde`

or
`hyperbolic`

.

`R`

— Dirichlet vectorsparse vector

Dirichlet vector, returned as a sparse vector. See Algorithms.

Typically, you use `R`

in a subsequent call to a solver
such as `assempde`

or
`hyperbolic`

.

This form of the finite element matrices eliminates Dirichlet
conditions from the problem using a linear algebra approach. The finite element
matrices reduce to the solution `u = B*(Kc\Fc) + ud`

, where
`B`

spans the null space of the columns of
`H`

(the Dirichlet condition matrix representing
*hu* = *r*). `R`

is the
Dirichlet condition vector for `Hu = R`

. `ud`

is the vector of boundary condition solutions for the Dirichlet conditions.
`u1`

= `Kc\Fc`

returns the
solution on the non-Dirichlet points.

See Systems of PDEs for details on the approach used to eliminate Dirichlet conditions.

This form of the finite element matrices converts Dirichlet
boundary conditions to Neumann boundary conditions using a stiff-spring
approximation. Using this approximation, `assempde`

returns a
matrix `Ks`

and a vector `Fs`

that represent
the combined finite element matrices. The approximate solution
`u`

is `u = Ks\Fs`

.

See Elliptic Equations. For details of the stiff-spring approximation, see Systems of PDEs.

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 Modeler app or through MATLAB

^{®}files. For 3-D geometry, import the geometry in STL file format.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

$$\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. | (1) |

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

`assema`

, while `assemb`

. When it is not necessary to distinguish
`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}*(

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

$$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

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

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}_{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

$$\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

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

$${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

To summarize, `assempde`

performs the following steps to obtain a
solution `u`

to an elliptic PDE:

Generate the finite element matrices [

`K`

,`M`

,`F`

,`Q`

,`G`

,`H`

,`R`

]. This step is equivalent to calling`assema`

to generate the matrices`K`

,`M`

, and`F`

, and also calling`assemb`

to generate the matrices`Q`

,`G`

,`H`

, and`R`

.Generate the combined finite element matrices [

`Kc`

,`Fc`

,`B`

,`ud`

]. The combined stiffness matrix is for the reduced linear system,`Kc = K + M + Q`

. The corresponding combined load vector is`Fc = F + G`

. The`B`

matrix spans the null space of the columns of`H`

(the Dirichlet condition matrix representing*hu*=*r*). The`R`

vector represents the Dirichlet conditions in`Hu = R`

. The`ud`

vector represents boundary condition solutions for the Dirichlet conditions.Calculate the solution

`u`

via`u = B*(Kc\Fc) + ud`

.

`assempde`

uses one of two algorithms for assembling a problem into
combined finite element matrix form. A *reduced linear system* form
leads to immediate solution via linear algebra. You choose the algorithm by the number of
outputs. For the reduced linear system form, request four outputs:

`[Kc,Fc,B,ud] = assempde(_)`

For the *stiff-spring approximation*, request two outputs:

`[Ks,Fs] = assempde(_)`

For details, see Reduced Linear System and Stiff-Spring Approximation.

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-*D*-by-*D*
tensor, and *D* is the geometry dimensions, 2 or 3.

For 2-D systems, the notation $$\nabla \cdot (c\otimes \nabla u)$$ represents
an *N*-by-1 matrix with an (*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$$

For 3-D systems, the notation $$\nabla \cdot (c\otimes \nabla u)$$ represents
an *N*-by-1 matrix with an (*i*,1)-component

$$\begin{array}{l}{\displaystyle \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 x}{c}_{i,j,1,3}\frac{\partial}{\partial z}\right){u}_{j}}\\ +{\displaystyle \sum _{j=1}^{N}\left(\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}+\frac{\partial}{\partial y}{c}_{i,j,2,3}\frac{\partial}{\partial z}\right){u}_{j}}\\ +{\displaystyle \sum _{j=1}^{N}\left(\frac{\partial}{\partial z}{c}_{i,j,3,1}\frac{\partial}{\partial x}+\frac{\partial}{\partial z}{c}_{i,j,3,2}\frac{\partial}{\partial y}+\frac{\partial}{\partial z}{c}_{i,j,3,3}\frac{\partial}{\partial z}\right){u}_{j}}\end{array}$$

The symbols **a** and **d** denote *N*-by-*N* matrices,
and **f** denotes a column vector of
length *N*.

The elements *c _{ijkl}*,

`c`

, `a`

, `d`

, and
`f`

. The case of identity, diagonal, and symmetric matrices are handled
as special cases. For the tensor 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}$$

For 2-D systems, the notation $$n\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\left(c\otimes \nabla u\right)$$ represents an *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)$$.

For 3-D systems, the notation $$n\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\left(c\otimes \nabla u\right)$$ represents an *N*-by-1 matrix with
(*i*,1)-component

$$\begin{array}{l}{\displaystyle \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{cos}(\alpha ){c}_{i,j,1,3}\frac{\partial}{\partial z}\right)}{u}_{j}\\ +{\displaystyle \sum _{j=1}^{N}\left(\mathrm{cos}(\beta ){c}_{i,j,2,1}\frac{\partial}{\partial x}+\mathrm{cos}(\beta ){c}_{i,j,2,2}\frac{\partial}{\partial y}+\mathrm{cos}(\beta ){c}_{i,j,2,3}\frac{\partial}{\partial z}\right)}{u}_{j}\\ +{\displaystyle \sum _{j=1}^{N}\left(\mathrm{cos}(\gamma ){c}_{i,j,3,1}\frac{\partial}{\partial x}+\mathrm{cos}(\gamma ){c}_{i,j,3,2}\frac{\partial}{\partial y}+\mathrm{cos}(\gamma ){c}_{i,j,3,3}\frac{\partial}{\partial z}\right)}{u}_{j}\end{array}$$

where the outward normal to the boundary is

$$n=\left(\mathrm{cos}\left(\alpha \right),\mathrm{cos}\left(\beta \right),\mathrm{cos}\left(\gamma \right)\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 mesh. Then the number of unknowns is

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

*B*´ *KBV* = *B*´
*F* –
*B*´*Ku _{d}*

which is symmetric and positive definite if *K* is.

The finite element method for 3-D geometry is similar to the 2-D method described in Elliptic Equations. The main difference is that the elements in 3-D geometry are tetrahedra, which means that the basis functions are different from those in 2-D geometry.

It is convenient to map a tetrahedron to a canonical tetrahedron with a local coordinate system (*r*,*s*,*t*).

In local coordinates, the point *p*1 is at (0,0,0), *p*2 is at (1,0,0), *p*3 is at (0,1,0), and *p*4 is at (0,0,1).

For a linear tetrahedron, the basis functions are

$$\begin{array}{c}{\varphi}_{1}=1-r-s-t\\ {\varphi}_{2}=r\\ {\varphi}_{3}=s\\ {\varphi}_{4}=t\end{array}$$

For a quadratic tetrahedron, there are additional nodes at the edge midpoints.

The corresponding basis functions are

$$\begin{array}{c}{\varphi}_{1}=2{\left(1-r-s-t\right)}^{2}-\left(1-r-s-t\right)\\ {\varphi}_{2}=2{r}^{2}-r\\ {\varphi}_{3}=2{s}^{2}-s\\ {\varphi}_{4}=2{t}^{2}-t\\ {\varphi}_{5}=4r\left(1-r-s-t\right)\\ {\varphi}_{6}=4rs\\ {\varphi}_{7}=4s\left(1-r-s-t\right)\\ {\varphi}_{8}=4t\left(1-r-s-t\right)\\ {\varphi}_{9}=4rt\\ {\varphi}_{10}=4st\end{array}$$

As in the 2-D case, a 3-D basis function *ϕ _{i}* takes the value 0 at all nodes

A modified version of this example exists on your system. Do you want to open this version instead?

You clicked a link that corresponds to this MATLAB command:

Run the command by entering it in the MATLAB Command Window. Web browsers do not support MATLAB commands.

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

Select web siteYou can also select a web site from the following list:

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

- América Latina (Español)
- Canada (English)
- United States (English)

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)