m, d, or a Coefficient for specifyCoefficients

Coefficients m, d, or a

This section describes how to write the m, d, or a coefficients in the system of equations

$$m\frac{{\partial }^{2}u}{\partial t^2}+d\frac{\partial u}{\partial t}-\nabla ·\left(c\otimes \nabla u\right)+au=f$$

or in the eigenvalue system

$$\begin{array}{l}-\nabla ·\left(c\otimes \nabla u\right)+au=\lambda du\\ \text{or}\\ -\nabla ·\left(c\otimes \nabla u\right)+au={\lambda }^{2}mu\end{array}$$

The topic applies to the recommended workflow for including coefficients in your model using specifyCoefficients.

If there are N equations in the system, then these coefficients represent N-by-N matrices.

For constant (numeric) coefficient matrices, represent each coefficient using a column vector with N² components. This column vector represents, for example, m(:).

For nonconstant coefficient matrices, see Nonconstant m, d, or a.

Note

The d coefficient takes a special matrix form when m is nonzero. See d Coefficient When m Is Nonzero.

Short m, d, or a vectors

Sometimes, your m, d, or a matrices are diagonal or symmetric. In these cases, you can represent m, d, or a using a smaller vector than one with N² components. The following sections give the possibilities.

• Scalar m, d, or a

• N-Element Column Vector m, d, or a

• N(N+1)/2-Element Column Vector m, d, or a

• N²-Element Column Vector m, d, or a

Scalar m, d, or a

The software interprets a scalar m, d, or a as a diagonal matrix.

$$\left(\begin{array}{cccc}a& 0& \cdots & 0\\ 0& a& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & a\end{array}\right)$$

N-Element Column Vector m, d, or a

The software interprets an N-element column vector m, d, or a as a diagonal matrix.

$$\left(\begin{array}{cccc}d(1)& 0& \cdots & 0\\ 0& d(2)& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & d(N)\end{array}\right)$$

N(N+1)/2-Element Column Vector m, d, or a

The software interprets an N(N+1)/2-element column vector m, d, or a as a symmetric matrix. In the following diagram, • means the entry is symmetric.

$$\left(\begin{array}{ccccc}a(1)& a(2)& a(4)& \cdots & a(N(N-1)/2)\\ •& a(3)& a(5)& \cdots & a(N(N-1)/2+1)\\ •& •& a(6)& \cdots & a(N(N-1)/2+2)\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ •& •& •& \cdots & a(N(N+1)/2)\end{array}\right)$$

Coefficient a(i,j) is in row (j(j–1)/2+i) of the vector a.

N²-Element Column Vector m, d, or a

The software interprets an N²-element column vector m, d, or a as a matrix.

$$\left(\begin{array}{cccc}d(1)& d(N+1)& \cdots & d(N^2-N+1)\\ d(2)& d(N+2)& \cdots & d(N^2-N+2)\\ ⋮& ⋮& \ddots & ⋮\\ d(N)& d(2N)& \cdots & d(N^2)\end{array}\right)$$

Coefficient a(i,j) is in row (N(j–1)+i) of the vector a.

Nonconstant m, d, or a

Note

If both m and d are nonzero, then d must be a constant scalar or vector, not a function.

If any of the m, d, or a coefficients is not constant, represent it as a function of the form

dcoeff = dcoeffunction(location,state)

Pass the coefficient to specifyCoefficients as a function handle, such as

specifyCoefficients(model,"d",@dcoeffunction,...)

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

dcoeff = ...
@(location,state) myfunWithAdditionalArgs(location,state,arg1,arg2...)
specifyCoefficients(model,d=dcoeff,...

• location is a structure with these fields:

  • location.x

  • location.y

  • location.z

  • location.subdomain

  The fields x, y, and z represent the x-, y-, and z- coordinates of points for which your function calculates coefficient values. The subdomain field represents the subdomain numbers, which currently apply only to 2-D models. The location fields are row vectors.

• state is a structure with these fields:

  • state.u

  • state.ux

  • state.uy

  • state.uz

  • state.time

  The state.u field represents the current value of the solution u. The state.ux, state.uy, and state.uz fields are estimates of the solution’s partial derivatives (∂u/∂x, ∂u/∂y, and ∂u/∂z) at the corresponding points of the location structure. The solution and gradient estimates are N-by-Nr matrices. The state.time field is a scalar representing time for time-dependent models.

Your function must return a matrix of size N1-by-Nr, where:

• N1 is the length of the vector representing the coefficient. There are several possible values of N1, detailed in Short m, d, or a vectors. 1 ≤ N1 ≤ N².

• Nr is the number of points in the location that the solver passes. Nr is equal to the length of the location.x or any other location field. The function should evaluate m, d, or a at these points.

For example, suppose N = 3, and you have 2-D geometry. Suppose your d matrix is of the form

$$d=\left[\begin{array}{ccc}1& s_1(x,y)& \sqrt{x^2+y^2}\\ s_1(x,y)& 4& -1\\ \sqrt{x^2+y^2}& -1& 9\end{array}\right]$$

where s₁(x,y) is 5 in subdomain 1, and is 10 in subdomain 2.

This d is a symmetric matrix. So it is natural to represent d as a N(N+1)/2-Element Column Vector m, d, or a:

$$\left(\begin{array}{ccccc}a(1)& a(2)& a(4)& \cdots & a(N(N-1)/2)\\ •& a(3)& a(5)& \cdots & a(N(N-1)/2+1)\\ •& •& a(6)& \cdots & a(N(N-1)/2+2)\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ •& •& •& \cdots & a(N(N+1)/2)\end{array}\right)$$

For that form, the following function is appropriate.

function dmatrix = dcoeffunction(location,state)

n1 = 6;
nr = numel(location.x);
dmatrix = zeros(n1,nr);
dmatrix(1,:) = ones(1,nr);
dmatrix(2,:) = 5*location.subdomain;
dmatrix(3,:) = 4*ones(1,nr);
dmatrix(4,:) = sqrt(location.x.^2 + location.y.^2);
dmatrix(5,:) = -ones(1,nr);
dmatrix(6,:) = 9*ones(1,nr);

To include this function as your d coefficient, pass the function handle @dcoeffunction:

specifyCoefficients(model,d=@dcoeffunction,...