## Specify Boundary Conditions

Before you create boundary conditions, you need to create a `PDEModel`

container. `PDEModel`

can accommodate one equation or a system of
*N* equations. For details, see Solve Problems Using PDEModel Objects.
Suppose that you have a container named `model`

, and that the geometry
is stored in `model`

. Examine the geometry to see the label of each
edge or face.

pdegplot(model,"EdgeLabels","on") % for 2-D pdegplot(model,"FaceLabels","on") % for 3-D

Now you can specify the boundary conditions for each edge or face. If you have a system of PDEs, you can set a different boundary condition for each component on each boundary edge or face.

If you do not specify a boundary condition for an edge or face, the default is the Neumann
boundary condition with the zero values for `"g"`

and
`"q"`

.

If the boundary condition is a function of position, time, or the solution
*u*, set boundary conditions by using the syntax in Nonconstant Boundary Conditions.

### Dirichlet Boundary Conditions

#### Scalar PDEs

The Dirichlet boundary condition implies that the solution *u* on
a particular edge or face satisfies the equation

*hu* = *r*,

where *h* and *r* are functions
defined on ∂Ω, and can be functions of space (*x*, *y*,
and, in 3-D, *z*), the solution *u*,
and, for time-dependent equations, time. Often, you take *h* = 1, and set *r* to
the appropriate value. You can specify Dirichlet boundary conditions
as the value of the solution `u`

on the boundary
or as a pair of the parameters `h`

and `r`

.

Suppose that you have a PDE model named `model`

,
and edges or faces `[e1,e2,e3]`

, where the solution *u* must
equal `2`

. Specify this boundary condition as follows.

% For 3-D geometry: applyBoundaryCondition(model,"dirichlet","Face",[e1,e2,e3],"u",2); % For 2-D geometry: applyBoundaryCondition(model,"dirichlet","Edge",[e1,e2,e3],"u",2);

If the solution on edges or faces `[e1,e2,e3]`

satisfies
the equation 2*u* = 3, specify the boundary condition
as follows.

% For 3-D geometry: applyBoundaryCondition(model,"dirichlet","Face",[e1,e2,e3],"r",3,"h",2); % For 2-D geometry: applyBoundaryCondition(model,"dirichlet","Edge",[e1,e2,e3],"r",3,"h",2);

If you do not specify

`"r"`

,`applyBoundaryCondition`

sets its value to`0`

.If you do not specify

`"h"`

,`applyBoundaryCondition`

sets its value to`1`

.

#### Systems of PDEs

The Dirichlet boundary condition for a system of PDEs is **hu** = **r**, where **h** is
a matrix, **u** is the solution vector,
and **r** is a vector.

Suppose that you have a PDE model named `model`

,
and edge or face labels `[e1,e2,e3]`

where the first
component of the solution *u* must equal `1`

,
while the second and third components must equal `2`

.
Specify this boundary condition as follows.

% For 3-D geometry: applyBoundaryCondition(model,"dirichlet","Face",[e1,e2,e3],... "u",[1,2,2],"EquationIndex",[1,2,3]); % For 2-D geometry: applyBoundaryCondition(model,"dirichlet","Edge",[e1,e2,e3],... "u",[1,2,2],"EquationIndex",[1,2,3]);

The

`"u"`

and`"EquationIndex"`

arguments must have the same length.If you exclude the

`"EquationIndex"`

argument, the`"u"`

argument must have length*N*.If you exclude the

`"u"`

argument,`applyBoundaryCondition`

sets the components in`"EquationIndex"`

to`0`

.

Suppose that you have a PDE model named `model`

,
and edge or face labels `[e1,e2,e3]`

where the first,
second, and third components of the solution *u* must
satisfy the equations 2*u _{1}* = 3, 4

*u*= 5, and 6

_{2}*u*= 7, respectively. Specify this boundary condition as follows.

_{3}H0 = [2 0 0; 0 4 0; 0 0 6]; R0 = [3;5;7]; % For 3-D geometry: applyBoundaryCondition(model,"dirichlet", ... "Face",[e1,e2,e3], ... "h",H0,"r",R0); % For 2-D geometry: applyBoundaryCondition(model,"dirichlet", ... "Edge",[e1,e2,e3], ... "h",H0,"r",R0);

The

`"r"`

parameter must be a numeric vector of length*N*. If you do not specify`"r"`

,`applyBoundaryCondition`

sets the values to`0`

.The

`"h"`

parameter can be an*N*-by-*N*numeric matrix or a vector of length*N*^{2}corresponding to the linear indexing form of the*N*-by-*N*matrix. For details about the linear indexing form, see Array Indexing. If you do not specify`"h"`

,`applyBoundaryCondition`

sets the value to the identity matrix.

### Neumann Boundary Conditions

#### Scalar PDEs

Generalized Neumann boundary conditions imply that the solution *u* on
the edge or face satisfies the equation

$$\overrightarrow{n}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\left(c\nabla u\right)+qu=g$$

The coefficient *c* is the same as the coefficient
of the second-order differential operator in the PDE equation

$$-\nabla \cdot \left(c\nabla u\right)+au=f\text{ondomain}\Omega $$

$$\overrightarrow{n}$$ is the outward unit normal. *q* and *g* are
functions defined on ∂Ω, and can be functions of space
(*x*, *y*, and, in 3-D, *z*),
the solution *u*, and, for time-dependent equations,
time.

Suppose that you have a PDE model named `model`

,
and edges or faces `[e1,e2,e3]`

where the solution *u* must
satisfy the Neumann boundary condition with `q = 2`

and ```
g
= 3
```

. Specify this boundary condition as follows.

% For 3-D geometry: applyBoundaryCondition(model,"neumann","Face",[e1,e2,e3],"q",2,"g",3); % For 2-D geometry: applyBoundaryCondition(model,"neumann","Edge",[e1,e2,e3],"q",2,"g",3);

If you do not specify

`"g"`

,`applyBoundaryCondition`

sets its value to`0`

.If you do not specify

`"q"`

,`applyBoundaryCondition`

sets its value to`0`

.

#### Systems of PDEs

Neumann boundary conditions for a system of PDEs is $$n\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\left(c\otimes \nabla u\right)+qu=g$$. For 2-D systems, the notation $$n\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\left(c\otimes \nabla u\right)$$ means the *N*-by-1
vector 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)\text{\hspace{0.17em}}}{u}_{j$$

where the outward normal vector of the boundary $$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)$$ means
the *N*-by-1 vector with (*i*,1)-component

$$\begin{array}{l}{\displaystyle \sum _{j=1}^{N}\left(\mathrm{sin}\left(\phi \right)\mathrm{cos}\left(\theta \right){c}_{i,j,1,1}\frac{\partial}{\partial x}+\mathrm{sin}\left(\phi \right)\mathrm{cos}\left(\theta \right){c}_{i,j,1,2}\frac{\partial}{\partial y}+\mathrm{sin}\left(\phi \right)\mathrm{cos}\left(\theta \right){c}_{i,j,1,3}\frac{\partial}{\partial z}\right){u}_{j}}\\ +{\displaystyle \sum _{j=1}^{N}\left(\mathrm{sin}\left(\phi \right)\mathrm{sin}\left(\theta \right){c}_{i,j,2,1}\frac{\partial}{\partial x}+\mathrm{sin}\left(\phi \right)\mathrm{sin}\left(\theta \right){c}_{i,j,2,2}\frac{\partial}{\partial y}+\mathrm{sin}\left(\phi \right)\mathrm{sin}\left(\theta \right){c}_{i,j,2,3}\frac{\partial}{\partial z}\right){u}_{j}}\\ +{\displaystyle \sum _{j=1}^{N}\left(\mathrm{cos}\left(\theta \right){c}_{i,j,3,1}\frac{\partial}{\partial x}+\mathrm{cos}\left(\theta \right){c}_{i,j,3,2}\frac{\partial}{\partial y}+\mathrm{cos}\left(\theta \right){c}_{i,j,3,3}\frac{\partial}{\partial z}\right){u}_{j}}\end{array}$$

where the outward normal vector of the boundary $$n=\left(\mathrm{sin}(\phi )\mathrm{cos}(\theta ),\mathrm{sin}(\phi )\mathrm{sin}(\theta ),\mathrm{cos}(\phi )\right)$$. For each edge or face segment,
there are a total of *N* boundary conditions.

Suppose that you have a PDE model named `model`

,
and edges or faces `[e1,e2,e3]`

where the first component
of the solution *u* must satisfy the Neumann boundary
condition with `q = 2`

and `g = 3`

,
and the second component must satisfy the Neumann boundary condition
with `q = 4`

and `g = 5`

. Specify
this boundary condition as follows.

Q = [2 0; 0 4]; G = [3;5]; % For 3-D geometry: applyBoundaryCondition(model,"neumann","Face",[e1,e2,e3],"q",Q,"g",G); % For 2-D geometry: applyBoundaryCondition(model,"neumann","Edge",[e1,e2,e3],"q",Q,"g",G);

The

`"g"`

parameter must be a numeric vector of length*N*. If you do not specify`"g"`

,`applyBoundaryCondition`

sets the values to`0`

.The

`"q"`

parameter can be an*N*-by-*N*numeric matrix or a vector of length*N*^{2}corresponding to the linear indexing form of the*N*-by-*N*matrix. For details about the linear indexing form, see Array Indexing. If you do not specify`"q"`

,`applyBoundaryCondition`

sets the values to`0`

.

### Mixed Boundary Conditions

If some equations in your system of PDEs must satisfy the Dirichlet boundary condition and
some must satisfy the Neumann boundary condition for the same geometric region,
use the `"mixed"`

parameter to apply boundary conditions in one
call. Note that `applyBoundaryCondition`

uses the default
Neumann boundary condition with `g = 0`

and ```
q =
0
```

for equations for which you do not explicitly specify a boundary
condition.

Suppose that you have a PDE model named `model`

,
and edge or face labels `[e1,e2,e3]`

where the first
component of the solution *u* must equal `11`

,
the second component must equal `22`

, and the third
component must satisfy the Neumann boundary condition with ```
q
= 3
```

and `g = 4`

. Express this boundary
condition as follows.

Q = [0 0 0; 0 0 0; 0 0 3]; G = [0;0;4]; % For 3-D geometry: applyBoundaryCondition(model,"mixed","Face",[e1,e2,e3],... "u",[11,22],"EquationIndex",[1,2],... "q",Q,"g",G); % For 2-D geometry: applyBoundaryCondition(model,"mixed",... "Edge",[e1,e2,e3],"u",[11,22],... "EquationIndex",[1,2],"q",Q,"g",G);

Suppose that you have a PDE model named `model`

,
and edges or faces `[e1,e2,e3]`

where the first component
of the solution *u* must satisfy the Dirichlet boundary
condition 2*u _{1}* = 3, the second
component must satisfy the Neumann boundary condition with

```
q
= 4
```

and `g = 5`

, and the third component
must satisfy the Neumann boundary condition with `q = 6`

and ```
g
= 7
```

. Express this boundary condition as follows.h = [2 0 0; 0 0 0; 0 0 0]; r = [3;0;0]; Q = [0 0 0; 0 4 0; 0 0 6]; G = [0;5;7]; % For 3-D geometry: applyBoundaryCondition(model,"mixed", ... "Face",[e1,e2,e3], ... "h",h,"r",r,"q",Q,"g",G); % For 2-D geometry: applyBoundaryCondition(model,"mixed", ... "Edge",[e1,e2,e3], ... "h",h,"r",r,"q",Q,"g",G);

### Nonconstant Boundary Conditions

Use functions to express nonconstant boundary conditions.

applyBoundaryCondition(model,"dirichlet", ... "Edge",1, ... "r",@myrfun); applyBoundaryCondition(model,"neumann", ... "Face",2, ... "g",@mygfun,"q",@myqfun); applyBoundaryCondition(model,"mixed", ... "Edge",[3,4], ... "u",@myufun, ... "EquationIndex",[2,3]);

Each function must have the following syntax.

`function bcMatrix = myfun(location,state)`

`solvepde`

or `solvepdeeig`

compute and
populate the data in the `location`

and `state`

structure arrays and pass this data to your function. You can define your function
so that its output depends on this data. You can use any names instead of
`location`

and `state`

, but the function
must have exactly two arguments.

`location`

— A structure containing the following fields. If you pass a name-value pair to`applyBoundaryCondition`

with`Vectorized`

set to`"on"`

, then`location`

can contain several evaluation points. If you do not set`Vectorized`

or use`Vectorized,"off"`

, then solvers pass just one evaluation point in each call.`location.x`

— The*x*-coordinate of the point or points`location.y`

— The*y*-coordinate of the point or points`location.z`

— For 3-D geometry, the*z*-coordinate of the point or points

Furthermore, if there are Neumann conditions, then solvers pass the following data in the

`location`

structure.`location.nx`

—*x*-component of the normal vector at the evaluation point or points`location.ny`

—*y*-component of the normal vector at the evaluation point or points`location.nz`

— For 3-D geometry,*z*-component of the normal vector at the evaluation point or points

`state`

— For transient or nonlinear problems.`state.u`

contains the solution vector at evaluation points.`state.u`

is an*N*-by-*M*matrix, where each column corresponds to one evaluation point, and*M*is the number of evaluation points.`state.time`

contains the time at evaluation points.`state.time`

is a scalar.

Your function returns `bcMatrix`

. This matrix
has the following form, depending on the boundary condition type.

`"u"`

—*N1*-by-*M*matrix, where each column corresponds to one evaluation point, and*M*is the number of evaluation points.*N1*is the length of the`"EquationIndex"`

argument. If there is no`"EquationIndex"`

argument, then*N1*=*N*.`"r"`

or`"g"`

—*N*-by-*M*matrix, where each column corresponds to one evaluation point, and*M*is the number of evaluation points.`"h"`

or`"q"`

—*N*^{2}-by-*M*matrix, where each column corresponds to one evaluation point via linear indexing of the underlying*N*-by-*N*matrix, and*M*is the number of evaluation points. Alternatively, an*N*-by-*N*-by-*M*array, where each evaluation point is an*N*-by-*N*matrix. For details about linear indexing, see Array Indexing.

If boundary conditions depend on `state.u`

or `state.time`

,
ensure that your function returns a matrix of `NaN`

of
the correct size when `state.u`

or `state.time`

are `NaN`

.
Solvers check whether a problem is nonlinear or time-dependent by
passing `NaN`

state values, and looking for returned `NaN`

values.

### Additional Arguments in Functions for Nonconstant Boundary Conditions

To use additional arguments in your function, wrap your function (that takes
additional arguments) with an anonymous function that takes only the
`location`

and `state`

arguments. For
example:

uVal = ... @(location,state) myfunWithAdditionalArgs(location,state,arg1,arg2...) applyBoundaryCondition(model,"mixed", ... "Edge",[3,4], ... "u",uVal, ... "EquationIndex",[2,3]);