Documentation

pdenonlin

Solve nonlinear elliptic PDE problem

pdenonlin IS NOT RECOMMENDED. Use solvepde instead.

Syntax

  • u = pdenonlin(model,c,a,f)
    example
  • u = pdenonlin(b,p,e,t,c,a,f)
    example
  • u = pdenonlin(___,Name,Value)
    example
  • [u,res] = pdenonlin(___)

Description

example

u = pdenonlin(model,c,a,f) solves the nonlinear PDE

(cu)+au=f,

with geometry, boundary conditions, and finite element mesh in model, and coefficients c, a, and f. In this context, "nonlinear" means some coefficient in c, a, or f depends on the solution u or its gradient. If the PDE is a system of equations (model.PDESystemSize > 1), then pdenonlin solves the system of equations

(cu)+au=f.

example

u = pdenonlin(b,p,e,t,c,a,f) solves the PDE with boundary conditions b, and finite element mesh (p,e,t).

example

u = pdenonlin(___,Name,Value), for any previous arguments, modifies the solution process with Name, Value pairs.

[u,res] = pdenonlin(___) also returns the norm of the Newton step residuals res.

Examples

Minimal Surface Problem

Solve a minimal surface problem. Because this problem has a nonlinear c coefficient, use pdenonlin to solve it.

Create a model and include circular geometry using the built-in circleg function.

model = createpde;
geometryFromEdges(model,@circleg);

Set the coefficients.

a = 0;
f = 0;
c = '1./sqrt(1+ux.^2+uy.^2)';

Set a Dirichlet boundary condition with value x2.

boundaryfun = @(region,state)region.x.^2;
applyBoundaryCondition(model,'edge',1:model.Geometry.NumEdges,...
    'u',boundaryfun,'Vectorized','on');

Generate a mesh and solve the problem.

generateMesh(model,'Hmax',0.1);
u = pdenonlin(model,c,a,f);
pdeplot(model,'xydata',u,'zdata',u);

Minimal Surface Problem Using [p,e,t] Mesh

Solve the minimal surface problem using the legacy approach for creating boundary conditions and geometry.

Create the geometry using the built-in circleg function. Plot the geometry to see the edge labels.

g = @circleg;
pdegplot(g,'EdgeLabels','on')
axis equal

Create Dirichlet boundary conditions with value x2. Create the following file and save it on your MATLAB® path. For details of this approach, see Boundary Conditions by Writing Functions.

function [qmatrix,gmatrix,hmatrix,rmatrix] = pdex2bound(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}
            hmatrix(k) = 1;
            hmatrix(k+ne) = 1;
            rmatrix(k) = x1^2;
            rmatrix(k+ne) = x2^2;
    end
end

Set the coefficients and boundary conditions.

a = 0;
f = 0;
c = '1./sqrt(1+ux.^2+uy.^2)';
b = @pdex2bound;

Generate a mesh and solve the problem.

[p,e,t] = initmesh(g,'Hmax',0.1);
u = pdenonlin(b,p,e,t,c,a,f);
pdeplot(p,e,t,'xydata',u,'zdata',u);

Nonlinear Problem with 3-D Geometry

Solve a nonlinear 3-D problem with nontrivial geometry.

Import the geometry from the BracketWithHole.stl file. Plot the geometry and face labels.

model = createpde();
importGeometry(model,'BracketWithHole.stl');
h = pdegplot(model,'FaceLabels','on');
h(1).FaceAlpha = 0.5;

Set a Dirichlet boundary condition with value 1000 on the back face, which is face 3. Set the large faces 1 and 5, and also the circular face 9, to have Neumann boundary conditions with value g = –10. Do not set boundary conditions on the other faces. Those faces default to Neumann boundary conditions with value g = 0.

applyBoundaryCondition(model,'face',3,'u',1000);
applyBoundaryCondition(model,'face',[1,5,9],'g',-10);

Set the c coefficient to 1, f to 0.1, and a to the nonlinear value '0.1 + 0.001*u.^2'.

c = 1;
f = 0.1;
a = '0.1 + 0.001*u.^2';

Generate the mesh and solve the PDE. Start from the initial guess u0 = 1000, which matches the value you set on face 3. Turn on the Report option to observe the convergence during the solution.

generateMesh(model);
u = pdenonlin(model,c,a,f,'U0',1000,'Report','on');
Iteration     Residual     Step size  Jacobian: full
   0          3.5985e-01
   1          1.0182e-01   1.0000000
   2          3.0318e-02   1.0000000
   3          8.7499e-03   1.0000000
   4          1.9023e-03   1.0000000
   5          1.5581e-04   1.0000000
   6          1.2424e-06   1.0000000

Plot the solution on the geometry boundary.

pdeplot3D(model,'colormapdata',u);

Input Arguments

collapse all

model — PDE modelPDEModel object

PDE model, specified as a PDEModel object.

Example: model = createpde(1)

c — PDE coefficientscalar or matrix | character array | coefficient function

PDE coefficient, specified as a scalar or matrix, as a character array, or as a coefficient function. c represents the c coefficient in the scalar PDE

(cu)+au=f,

or in the system of PDEs

(cu)+au=f.

You can specifyc in various ways, detailed in c Coefficient for Systems. See also Scalar PDE Coefficients, Specify Scalar PDE Coefficients in String Form, Specify 2-D Scalar Coefficients in Function Form, Specify 3-D PDE Coefficients in Function Form, and Coefficients for Systems of PDEs.

Example: 'cosh(x+y.^2)'

Data Types: double | char | function_handle
Complex Number Support: Yes

a — PDE coefficientscalar or matrix | character array | coefficient function

PDE coefficient, specified as a scalar or matrix, as a character array, or as a coefficient function. a represents the a coefficient in the scalar PDE

(cu)+au=f,

or in the system of PDEs

(cu)+au=f.

You can specifya in various ways, detailed in a or d Coefficient for Systems. See also Scalar PDE Coefficients, Specify Scalar PDE Coefficients in String Form, Specify 2-D Scalar Coefficients in Function Form, Specify 3-D PDE Coefficients in Function Form, and Coefficients for Systems of PDEs.

Example: 2*eye(3)

Data Types: double | char | function_handle
Complex Number Support: Yes

f — PDE coefficientscalar or matrix | character array | coefficient function

PDE coefficient, specified as a scalar or matrix, as a character array, or as a coefficient function. f represents the f coefficient in the scalar PDE

(cu)+au=f,

or in the system of PDEs

(cu)+au=f.

You can specifyf in various ways, detailed in f Coefficient for Systems. See also Scalar PDE Coefficients, Specify Scalar PDE Coefficients in String Form, Specify 2-D Scalar Coefficients in Function Form, Specify 3-D PDE Coefficients in Function Form, and Coefficients for Systems of PDEs.

Example: char('sin(x)';'cos(y)';'tan(z)')

Data Types: double | char | function_handle
Complex Number Support: Yes

b — Boundary conditionsboundary matrix | boundary file

Boundary conditions, specified as a boundary matrix or boundary file. Pass a boundary file as a function handle or as a string naming the file.

For more information on boundary conditions, see Forms of Boundary Condition Specification.

Example: b = 'circleb1' or equivalently b = @circleb1

Data Types: double | char | function_handle

p — Mesh nodesoutput of initmesh | output of meshToPet

Mesh nodes, specified as the output of initmesh or meshToPet. For the structure of a p matrix, see Mesh Data for [p,e,t] Triples: 2-D and Mesh Data for [p,e,t] Triples: 3-D.

Example: [p,e,t] = initmesh(g)

Data Types: double

e — Mesh edgesoutput of initmesh | output of meshToPet

Mesh edges, specified as the output of initmesh or meshToPet. For the structure of e, see Mesh Data for [p,e,t] Triples: 2-D and Mesh Data for [p,e,t] Triples: 3-D.

Example: [p,e,t] = initmesh(g)

Data Types: double

t — Mesh elementsoutput of initmesh | output of meshToPet

Mesh elements, specified as the output of initmesh or meshToPet. Mesh elements are the triangles or tetrahedra that form the finite element mesh. For the structure of a t matrix, see Mesh Data for [p,e,t] Triples: 2-D and Mesh Data for [p,e,t] Triples: 3-D.

Example: [p,e,t] = initmesh(g)

Data Types: double

Name-Value Pair Arguments

Specify optional comma-separated pairs of Name,Value arguments. Name is the argument name and Value is the corresponding value. Name must appear inside single quotes (' '). You can specify several name and value pair arguments in any order as Name1,Value1,...,NameN,ValueN.

Example:

'Jacobian' — Approximation of Jacobian'full' (3-D default) | 'fixed' (2-D default) | 'lumped'

Approximation of Jacobian, specified as 'full', 'fixed', or 'lumped'.

  • 'full' means numerical evaluation of the full Jacobian based on the sparse version of the numjac function. 3-D geometry uses only 'full', any other specification yields an error.

  • 'fixed' specifies a fixed-point iteration matrix where the Jacobian is approximated by the stiffness matrix. This is the 2-D geometry default.

  • 'lumped' specifies a "lumped" approximation as described in Nonlinear Equations. This approximation is based on the numerical differentiation of the coefficients.

Example: u = pdenonlin(model,c,a,f,'Jacobian','full')

Data Types: char

'U0' — Initial solution guess0 (default) | scalar | string | vector

Initial solution guess, specified as a scalar, string, or vector. For details, see Solve PDEs with Initial Conditions.

  • A scalar specifies a constant initial condition for either a scalar or PDE system.

  • For scalar problems, write a string using the same syntax as Specify Scalar PDE Coefficients in String Form.

  • For systems of N equations, write a character array with N rows, where each row has the syntax of Specify Scalar PDE Coefficients in String Form.

  • For systems of N equations, and a mesh with Np nodes, give a column vector with N*Np components. The nodes are either model.Mesh.Nodes, or the p data from initmesh or meshToPet. See Mesh Data.

    The first Np elements contain the values of component 1, where the value of element k corresponds to node p(k). The next Np points contain the values of component 2, etc. It can be convenient to first represent the initial conditions u0 as an Np-by-N matrix, where the first column contains entries for component 1, the second column contains entries for component 2, etc. The final representation of the initial conditions is u0(:).

Example: u = pdenonlin(model,c,a,f,'U0','x.^2-y.^2')

Data Types: double | char
Complex Number Support: Yes

'Tol' — Residual size at termination1e-4 (default) | positive scalar

Residual size at termination, specified as a positive scalar. pdenonlin iterates until the residual size is less than 'Tol'.

Example: u = pdenonlin(model,c,a,f,'Tol',1e-6)

Data Types: double

'MaxIter' — Maximum number of Gauss-Newton iterations25 (default) | positive integer

Maximum number of Gauss-Newton iterations, specified as a positive integer.

Example: u = pdenonlin(model,c,a,f,'MaxIter',12)

Data Types: double

'MinStep' — Minimum damping of search direction1/2^16 (default) | positive scalar

Minimum damping of search direction, specified as a positive scalar.

Example: u = pdenonlin(model,c,a,f,'MinStep',1e-3)

Data Types: double

'Report' — Print convergence information'off' (default) | 'on'

Print convergence information, specified as 'off' or 'on'.

Example: u = pdenonlin(model,c,a,f,'Report','on')

Data Types: char

'Norm' — Residual normInf (default) | p value for Lp norm | 'energy'

Residual norm, specified as the p value for Lp norm, or as the string 'energy'. p can be any positive real value, Inf, or -Inf. The p norm of a vector v is sum(abs(v)^p)^(1/p). See norm.

Example: u = pdenonlin(model,c,a,f,'Norm',2)

Data Types: double | char

Output Arguments

collapse all

u — PDE solutionvector

PDE solution, returned as a vector.

  • If the PDE is scalar, meaning only one equation, then u is a column vector representing the solution u at each node in the mesh. u(i) is the solution at the ith column of model.Mesh.Nodes or the ith column of p.

  • If the PDE is a system of N > 1 equations, then u is a column vector with N*Np elements, where Np is the number of nodes in the mesh. The first Np elements of u represent the solution of equation 1, then next Np elements represent the solution of equation 2, etc.

To obtain the solution at an arbitrary point in the geometry, use pdeInterpolant.

To plot the solution, use pdeplot for 2-D geometry, or see Plot 3-D Solutions and Their Gradients.

res — Norm of Newton step residualsscalar

Norm of Newton step residuals, returned as a scalar. For information about the algorithm, see Nonlinear Equations.

More About

collapse all

Tips

  • If the Newton iteration does not converge, pdenonlin displays the error message Too many iterations or Stepsize too small.

  • If the initial guess produces matrices containing NaN or Inf elements, pdenonlin displays the error message Unsuitable initial guess U0 (default: U0 = 0).

  • If you have very small coefficients, or very small geometric dimensions, pdenonlin can fail to converge, or can converge to an incorrect solution. If so, you can sometimes obtain better results by scaling the coefficients or geometry dimensions to be of order one.

Algorithms

The pdenonlin algorithm is a Gauss-Newton iteration scheme applied to the finite element matrices. For details, see Nonlinear Equations.

See Also

Introduced before R2006a

Was this topic helpful?