Documentation |
This section describes how to write the coefficient f in the equation
$$-\nabla \cdot \left(c\otimes \nabla u\right)+au=f,$$
or in similar equations. The number of rows in f indicates N, the number of equations, see Systems of PDEs. Give f as any of the following:
A column vector with N components. For example, if N = 3, f could be:
f = [3;4;10];
A character array with N rows. The rows of the character array are MATLAB^{®} expressions as described in Scalar PDE Coefficients in String Form, with additional options for nonlinear equations. The additional options are:
Represent the ith component of the solution u using 'u(i)'.
Similarly, represent the ith components of the gradients of the solution u using 'ux(i)' and 'uy(i)'.
Pad the rows with spaces so each row has the same number of characters (char does this automatically). For example, if N = 3, f could be:
f = char('sin(x)+cos(y)','cosh(x.*y)*(1+u(1).^2)','x.*y./(1+x.^2+y.^2)')
f = sin(x)+cos(y) cosh(x.*y)*(1+u(1).^2) x.*y./(1+x.^2+y.^2)
A function of the form as described in Scalar PDE Coefficients in Function Form. The function should return a matrix of size N-by-Nt, where Nt is the number of triangles in the mesh. The function should evaluate f at the triangle centroids, as in Scalar PDE Coefficients in Function Form. Give solvers the function name as a string 'filename', or as a function handle @filename, where filename.m is a file on your MATLAB path. For details on writing the function, see Calculate Coefficients in Function Form.
For example, if N = 3, f could be:
function f = fcoeffunction(p,t,u,time) N = 3; % Number of equations % Triangle point indices it1 = t(1,:); it2 = t(2,:); it3 = t(3,:); % Find centroids of triangles xpts = (p(1,it1)+p(1,it2)+p(1,it3))/3; ypts = (p(2,it1)+p(2,it2)+p(2,it3))/3; [ux,uy] = pdegrad(p,t,u); % Approximate derivatives uintrp = pdeintrp(p,t,u); % Interpolated values at centroids nt = size(t,2); % Number of columns f = zeros(N,nt); % Allocate f % Now the particular functional form of f f(1,:) = xpts - ypts + uintrp(1,:); f(2,:) = 1 + tanh(ux(1,:)) + tanh(uy(3,:)); f(3,:) = (5+uintrp(3,:)).*sqrt(xpts.^2+ypts.^2);
Because this function depends on the solution u, if the equation is elliptic, use the pdenonlin solver. The initial value can be all 0s in the case of Dirichlet boundary conditions:
np = size(p,2); % number of points u0 = zeros(N*np,1); % initial guess
Caution It is not reliable to specify f as a scalar or single string expression. Sometimes the toolbox can expand the single input to a vector or character array with N identical rows. But you can get an error when the toolbox fails to determine N. Instead of a scalar or single string, for reliability specify f as a column vector or character array with N rows. |