**THIS PAGE DESCRIBES
THE LEGACY WORKFLOW.** New features might not be
compatible with the legacy workflow. For the corresponding step in
the recommended workflow, see Specify Boundary Conditions.

This section shows how to express boundary conditions for 2-D geometry using the legacy function syntax. However, the recommended way to express boundary conditions is to use Specify Boundary Conditions.

To use this legacy syntax, write the functions using the templates in Boundary Conditions for Scalar PDE or Boundary Conditions for PDE Systems.

For a scalar PDE, some boundary segments can have Dirichlet conditions, and some boundary segments can have generalized Neumann conditions.

Dirichlet boundary conditions are

*hu* = *r*,

where *h* and *r* can be functions
of *x*, *y*, the solution *u*,
the edge segment index, and, for parabolic and hyperbolic equations,
time.

Generalized Neumann boundary conditions are $$\overrightarrow{n}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\left(c\nabla u\right)+qu=g$$ on ∂Ω.

$$\overrightarrow{n}$$ is the outward unit normal. *g* and *q* are
functions defined on ∂Ω, and can be functions of *x*, *y*,
the solution *u*, the edge segment index, and, for
parabolic and hyperbolic equations, time.

To write a function file, say `pdebound.m`

,
use the following syntax:

[qmatrix,gmatrix,hmatrix,rmatrix] = pdebound(p,e,u,time)

Your function returns matrices `qmatrix`

, `gmatrix`

, `hmatrix`

,
and `rmatrix`

, based on these inputs:

If your boundary conditions do not depend on `u`

or `time`

,
those inputs are `[]`

. If your boundary conditions
do depend on `u`

or `time`

, then
when `u`

or `time`

are `NaN`

,
ensure that the outputs such as `qmatrix`

consist
of matrices of `NaN`

of the correct size. This signals
to solvers, such as `parabolic`

, to use a time-dependent
or solution-dependent algorithm.

Before specifying boundary conditions, you need to know the boundary labels. See Identify Boundary Labels.

The PDE solver, such as `assempde`

or `adaptmesh`

,
passes a matrix `p`

of points and `e`

of
edges. `e`

has seven rows and `ne`

columns,
where you do not necessarily know in advance the size `ne`

.

`p`

is a 2-by-*N*matrix, where_{p}`p(1,k)`

is the*x*-coordinate of point`k`

, and`p(2,k)`

is the*y*-coordinate of point`k`

.`e`

is a 7-by-`ne`

matrix, where`e(1,k)`

is the index of the first point of edge`k`

.`e(2,k)`

is the index of the second point of edge`k`

.`e(5,k)`

is the label of the geometry edge of edge`k`

(see Identify Boundary Labels).

`e`

contains an entry for every finite element edge that lies on an exterior boundary.

Use the following template for your boundary file.

function [qmatrix,gmatrix,hmatrix,rmatrix] = pdebound(p,e,u,time) ne = size(e,2); % number of edges qmatrix = zeros(1,ne); gmatrix = qmatrix; hmatrix = zeros(1,2*ne); rmatrix = hmatrix; for k = 1:ne x1 = p(1,e(1,k)); % x at first point in segment x2 = p(1,e(2,k)); % x at second point in segment xm = (x1 + x2)/2; % x at segment midpoint y1 = p(2,e(1,k)); % y at first point in segment y2 = p(2,e(2,k)); % y at second point in segment ym = (y1 + y2)/2; % y at segment midpoint switch e(5,k) case {some_edge_labels} % Fill in hmatrix,rmatrix or qmatrix,gmatrix case {another_list_of_edge_labels} % Fill in hmatrix,rmatrix or qmatrix,gmatrix otherwise % Fill in hmatrix,rmatrix or qmatrix,gmatrix end end

For each column `k`

in `e`

,
entry `k`

of `rmatrix`

is the value
of `rmatrix`

at the first point in the edge, and
entry `ne`

+ `k`

is
the value at the second point in the edge. For example, if *r* = *x*^{2} + *y*^{4},
then write these lines:

rmatrix(k) = x1^2 + y1^4; rmatrix(k+ne) = x2^2 + y2^4;

The syntax for `hmatrix`

is identical: entry `k`

of `hmatrix`

is
the value of `r`

at the first point in the edge,
and entry `k`

+ `ne`

is
the value at the second point in the edge.

For each column `k`

in `e`

,
entry `k`

of `qmatrix`

is the value
of `qmatrix`

at the midpoint in the edge. For example,
if *q* = *x*^{2} + *y*^{4},
then write these lines:

qmatrix(k) = xm^2 + ym^4;

The syntax for `gmatrix`

is identical: entry `k`

of `gmatrix`

is
the value of `gmatrix`

at the midpoint in the edge.

If the coefficients depend on the solution `u`

,
use the element `u(e(1,k))`

as the solution value
at the first point of edge `k`

, and `u(e(2,k))`

as
the solution value at the second point of edge `k`

.

For example, consider the following geometry, a rectangle with a circular hole.

% Rectangle is code 3, 4 sides, % followed by x-coordinates and then y-coordinates R1 = [3,4,-1,1,1,-1,-.4,-.4,.4,.4]'; % Circle is code 1, center (.5,0), radius .2 C1 = [1,.5,0,.2]'; % Pad C1 with zeros to enable concatenation with R1 C1 = [C1;zeros(length(R1)-length(C1),1)]; geom = [R1,C1]; % Names for the two geometric objects ns = (char('R1','C1'))'; % Set formula sf = 'R1-C1'; % Create geometry gd = decsg(geom,sf,ns); % View geometry pdegplot(gd,'EdgeLabels','on') xlim([-1.1 1.1]) axis equal

Suppose the boundary conditions on the outer boundary (segments
1 through 4) are Dirichlet, with the value *u*(*x*,*y*) = *t*(*x* – *y*),
where *t* is time. Suppose the circular boundary
(segments 5 through 8) has a generalized Neumann condition, with *q* = 1 and *g* = *x*^{2} + *y*^{2}.

Write the following boundary file to represent the boundary conditions:

function [qmatrix,gmatrix,hmatrix,rmatrix] = pdebound(p,e,u,time) ne = size(e,2); % number of edges qmatrix = zeros(1,ne); gmatrix = qmatrix; hmatrix = zeros(1,2*ne); rmatrix = hmatrix; for k = 1:ne x1 = p(1,e(1,k)); % x at first point in segment x2 = p(1,e(2,k)); % x at second point in segment xm = (x1 + x2)/2; % x at segment midpoint y1 = p(2,e(1,k)); % y at first point in segment y2 = p(2,e(2,k)); % y at second point in segment ym = (y1 + y2)/2; % y at segment midpoint switch e(5,k) case {1,2,3,4} % rectangle boundaries hmatrix(k) = 1; hmatrix(k+ne) = 1; rmatrix(k) = time*(x1 - y1); rmatrix(k+ne) = time*(x2 - y2); otherwise % same as case {5,6,7,8}, circle boundaries qmatrix(k) = 1; gmatrix(k) = xm^2 + ym^2; end end

The general mixed-boundary conditions for PDE systems of *N* equations
(see Equations You Can Solve Using Legacy Functions) are

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

The notation $$n\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\left(c\otimes \nabla u\right)$$ means the *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)\text{\hspace{0.17em}}}{u}_{j$$

where the outward normal vector of the boundary $$n=\left(\mathrm{cos}(\alpha ),\mathrm{sin}(\alpha )\right)$$. For each edge segment 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 solver computes Lagrange
multipliers *µ* such that the Dirichlet conditions
are satisfied.

To write a function file, say `pdebound.m`

,
use the following syntax:

[qmatrix,gmatrix,hmatrix,rmatrix] = pdebound(p,e,u,time)

Your function returns matrices `qmatrix`

, `gmatrix`

, `hmatrix`

,
and `rmatrix`

, based on these inputs:

If your boundary conditions do not depend on `u`

or `time`

,
those inputs are `[]`

. If your boundary conditions
do depend on `u`

or `time`

, then
when `u`

or `time`

are `NaN`

,
ensure that the outputs such as `qmatrix`

consist
of matrices of `NaN`

of the correct size. This signals
to solvers, such as `parabolic`

, to use a time-dependent
or solution-dependent algorithm.

Before specifying boundary conditions, you need to know the boundary labels. See Identify Boundary Labels.

A PDE solver, such as `assempde`

or `adaptmesh`

,
passes a matrix `p`

of points and `e`

of
edges. `e`

has seven rows and `ne`

columns,
where you do not necessarily know in advance the size `ne`

.

`p`

is a 2-by-*N*matrix, where_{p}`p(1,k)`

is the*x*-coordinate of point`k`

, and`p(2,k)`

is the*y*-coordinate of point`k`

.`e`

is a 7-by-`ne`

matrix, where`e(1,k)`

is the index of the first point of edge`k`

.`e(2,k)`

is the index of the second point of edge`k`

.`e(5,k)`

is the label of the geometry edge of edge`k`

(see Identify Boundary Labels).

`e`

contains an entry for every finite element edge that lies on an exterior boundary.

Let *N* be the dimension of the system of PDEs;
see Equations You Can Solve Using Legacy Functions. Use
the following template for your boundary file.

function [qmatrix,gmatrix,hmatrix,rmatrix] = pdebound(p,e,u,time) N = 3; % Set N = the number of equations ne = size(e,2); % number of edges qmatrix = zeros(N^2,ne); gmatrix = zeros(N,ne); hmatrix = zeros(N^2,2*ne); rmatrix = zeros(N,2*ne); for k = 1:ne x1 = p(1,e(1,k)); % x at first point in segment x2 = p(1,e(2,k)); % x at second point in segment xm = (x1 + x2)/2; % x at segment midpoint y1 = p(2,e(1,k)); % y at first point in segment y2 = p(2,e(2,k)); % y at second point in segment ym = (y1 + y2)/2; % y at segment midpoint switch e(5,k) case {some_edge_labels} % Fill in hmatrix,rmatrix or qmatrix,gmatrix case {another_list_of_edge_labels} % Fill in hmatrix,rmatrix or qmatrix,gmatrix otherwise % Fill in hmatrix,rmatrix or qmatrix,gmatrix end end

For the boundary file, you represent the matrix **h** for each edge segment as a vector, taking
the matrix column-wise, as `hmatrix(:)`

. Column `k`

of `hmatrix`

corresponds
to the matrix at the first edge point `e(1,k)`

, and
column `k`

+ `ne`

corresponds
to the matrix at the second edge point `e(2,k)`

.

Similarly, you represent each vector **r** for
an edge as a column in the matrix `rmatrix`

. Column `k`

corresponds
to the vector at the first edge point `e(1,k)`

, and
column `k`

+ `ne`

corresponds
to the vector at the second edge point `e(2,k)`

.

Represent the entries for the matrix **q** for
each edge segment as a vector, `qmatrix(:)`

, similar
to the matrix `hmatrix(:)`

. Similarly, represent **g** for each edge segment is a column vector
in the matrix `gmatrix`

. Unlike **h** and **r**, which have two columns for each segment, **q** and **g** have
just one column for each segment, which is the value of the function
at the midpoint of the edge segment.

For example, consider the following geometry, a rectangle with a circular hole.

% Rectangle is code 3, 4 sides, % followed by x-coordinates and then y-coordinates R1 = [3,4,-1,1,1,-1,-.4,-.4,.4,.4]'; % Circle is code 1, center (.5,0), radius .2 C1 = [1,.5,0,.2]'; % Pad C1 with zeros to enable concatenation with R1 C1 = [C1;zeros(length(R1)-length(C1),1)]; geom = [R1,C1]; % Names for the two geometric objects ns = (char('R1','C1'))'; % Set formula sf = 'R1-C1'; % Create geometry gd = decsg(geom,sf,ns); % View geometry pdegplot(gd,'EdgeLabels','on') xlim([-1.1 1.1]) axis equal

Suppose *N* = 3. Suppose the boundary conditions
are mixed. There is `M`

= 1 Dirichlet condition:

The first component of

*u*= 0 on the rectangular segments (numbers 1–4). So**h**(1,1) = 1 and**r**(1) = 0 for those segments.The second components of

*u*= 0 on the circular segments (numbers 5–8). So**h**(2,2) = 1 and**r**(2) = 0 for those segments.On the rectangular segments (numbers 1–4),

$$q=\left(\begin{array}{ccc}0& 1& 1\\ 0& 0& 0\\ 1& 1& 0\end{array}\right)$$

and

$$g=\left(\begin{array}{c}1+{x}^{2}\\ 0\\ 1+{y}^{2}\end{array}\right)$$

On the circular segments (numbers 5–8),

$$q=\left(\begin{array}{ccc}0& 1+{x}^{2}& 2+{y}^{2}\\ 0& 0& 0\\ 1+{x}^{4}& 1+{y}^{4}& 0\end{array}\right)$$

and

$$g=\left(\begin{array}{c}\mathrm{cos}(\pi x)\\ 0\\ \mathrm{tanh}(x+y)\end{array}\right)$$

Write the following boundary file to represent the boundary conditions:

function [qmatrix,gmatrix,hmatrix,rmatrix] = pdebound(p,e,u,time) N = 3; ne = size(e,2); % number of edges qmatrix = zeros(N^2,ne); gmatrix = zeros(N,ne); hmatrix = zeros(N^2,2*ne); rmatrix = zeros(N,2*ne); for k = 1:ne x1 = p(1,e(1,k)); % x at first point in segment x2 = p(1,e(2,k)); % x at second point in segment xm = (x1 + x2)/2; % x at segment midpoint y1 = p(2,e(1,k)); % y at first point in segment y2 = p(2,e(2,k)); % y at second point in segment ym = (y1 + y2)/2; % y at segment midpoint switch e(5,k) case {1,2,3,4} hk = zeros(N); hk(1,1) = 1; hk = hk(:); hmatrix(:,k) = hk; hmatrix(:,k+ne) = hk; rk = zeros(N,1); % Not strictly necessary rmatrix(:,k) = rk; % These are already 0 rmatrix(:,k+ne) = rk; qk = zeros(N); qk(1,2) = 1; qk(1,3) = 1; qk(3,1) = 1; qk(3,2) = 1; qk = qk(:); qmatrix(:,k) = qk; gk = zeros(N,1); gk(1) = 1+xm^2; gk(3) = 1+ym^2; gmatrix(:,k) = gk; case {5,6,7,8} hk = zeros(N); hk(2,2) = 1; hk = hk(:); hmatrix(:,k) = hk; hmatrix(:,k+ne) = hk; rk = zeros(N,1); % Not strictly necessary rmatrix(:,k) = rk; % These are already 0 rmatrix(:,k+ne) = rk; qk = zeros(N); qk(1,2) = 1+xm^2; qk(1,3) = 2+ym^2; qk(3,1) = 1+xm^4; qk(3,2) = 1+ym^4; qk = qk(:); qmatrix(:,k) = qk; gk = zeros(N,1); gk(1) = cos(pi*xm); gk(3) = tanh(xm*ym); gmatrix(:,k) = gk; end end