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 2u = 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 2u₁ = 3, 4u₂ = 5, and 6u₃ = 7, respectively. Specify this boundary condition as follows.

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² 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}·\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{ on domain }\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·\left(c\otimes \nabla u\right)+qu=g$$. For example, in case of circumferential and spherical boundaries, the generalized versions of the Neumann boundary condition are as follows:

• If the boundary is a circumference (2-D case), the outward normal vector of the boundary of the boundary is given by $$n=\left(\cos (\phi ),\sin (\phi )\right)$$, the notation $$n·\left(c\otimes \nabla u\right)$$ means the N-by-1 vector, for which the (i,1)-component is as follows:

  $${\displaystyle \sum _{j=1}^N\left(\cos (\phi )c_{i,j,1,1}\frac{\partial }{\partial x}+\cos (\phi )c_{i,j,1,2}\frac{\partial }{\partial y}+\sin (\phi )c_{i,j,2,1}\frac{\partial }{\partial x}+\sin (\phi )c_{i,j,2,2}\frac{\partial }{\partial y}\right)}{u}_j$$

  

• If the boundary is a spherical surface (3-D case), than the outward normal vector of the boundary is given by $$n=\left(\sin (\theta )\cos (\phi ),\sin (\theta )\sin (\phi ),\cos (\theta )\right)$$, and the notation $$n·\left(c\otimes \nabla u\right)$$ means the N-by-1 vector, for which the (i,1)-component is as follows:

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

  

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² 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 2u₁ = 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²-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.

See Specify Nonconstant Boundary Conditions.

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]);