Documentation

This is machine translation

Translated by Microsoft
Mouse over text to see original. Click the button below to return to the English verison of the page.

pdeeig

Solve eigenvalue PDE problem

pdeeig is not recommended. Use solvepdeeig instead.

Syntax

  • [v,l] = pdeeig(model,c,a,d,r)
    example
  • [v,l] = pdeeig(b,p,e,t,c,a,d,r)
    example
  • [v,l] = pdeeig(Kc,B,M,r)
    example

Description

example

[v,l] = pdeeig(model,c,a,d,r) produces the solution to the FEM formulation of the scalar PDE eigenvalue problem

(cu)+au=λdu on Ω

or the system PDE eigenvalue problem

(cu)+au=λdu on Ω

with geometry, boundary conditions, and mesh specified in model, a PDEModel object. See Solve Problems Using Legacy PDEModel Objects.

The eigenvalue PDE problem is a homogeneous problem, i.e., only boundary conditions where g = 0 and r = 0 can be used. The nonhomogeneous part is removed automatically.

example

[v,l] = pdeeig(b,p,e,t,c,a,d,r) solves for boundary conditions described in b, and the finite element mesh in [p,e,t].

example

[v,l] = pdeeig(Kc,B,M,r) produces the solution to the generalized sparse matrix eigenvalue problem

Kc ui = λB´MBui
u = Bui

with Real(λ) in the interval r.

Examples

collapse all

Compute the eigenvalues that are less than 100, and compute the corresponding eigenmodes for $- \nabla u = \lambda u$ on the geometry of the L-shaped membrane.

model = createpde;
geometryFromEdges(model,@lshapeg);
applyBoundaryCondition(model,'edge',1:model.Geometry.NumEdges,'u',0);
generateMesh(model,'Hmax',0.02);
c = 1;
a = 0;
d = 1;
r = [-Inf 100];
[v,l] = pdeeig(model,c,a,d,r);
l(1)                    % first eigenvalue
              Basis= 10,  Time=   3.18,  New conv eig=  0
              Basis= 11,  Time=   3.35,  New conv eig=  0
              Basis= 12,  Time=   3.58,  New conv eig=  0
              Basis= 13,  Time=   3.77,  New conv eig=  0
              Basis= 14,  Time=   3.98,  New conv eig=  0
              Basis= 15,  Time=   4.16,  New conv eig=  1
              Basis= 16,  Time=   4.29,  New conv eig=  1
              Basis= 17,  Time=   4.48,  New conv eig=  3
              Basis= 18,  Time=   4.69,  New conv eig=  4
              Basis= 19,  Time=   4.87,  New conv eig=  4
              Basis= 20,  Time=   5.05,  New conv eig=  4
              Basis= 21,  Time=   5.13,  New conv eig=  4
              Basis= 22,  Time=   5.32,  New conv eig=  4
              Basis= 23,  Time=   5.42,  New conv eig=  4
              Basis= 24,  Time=   5.60,  New conv eig=  5
              Basis= 25,  Time=   5.78,  New conv eig=  5
              Basis= 26,  Time=   5.99,  New conv eig=  5
              Basis= 27,  Time=   6.11,  New conv eig=  6
              Basis= 28,  Time=   6.30,  New conv eig=  7
              Basis= 29,  Time=   6.48,  New conv eig=  7
              Basis= 30,  Time=   6.59,  New conv eig=  7
              Basis= 31,  Time=   6.78,  New conv eig=  7
              Basis= 32,  Time=   6.99,  New conv eig=  8
              Basis= 33,  Time=   7.11,  New conv eig=  8
              Basis= 34,  Time=   7.42,  New conv eig=  8
              Basis= 35,  Time=   7.69,  New conv eig=  9
              Basis= 36,  Time=   7.84,  New conv eig=  9
              Basis= 37,  Time=   8.09,  New conv eig=  9
              Basis= 38,  Time=   8.37,  New conv eig=  9
              Basis= 39,  Time=   8.55,  New conv eig=  9
              Basis= 40,  Time=   8.82,  New conv eig=  9
              Basis= 41,  Time=   9.15,  New conv eig=  9
              Basis= 42,  Time=   9.46,  New conv eig= 10
              Basis= 43,  Time=   9.65,  New conv eig= 10
              Basis= 44,  Time=   9.96,  New conv eig= 12
              Basis= 45,  Time=  10.23,  New conv eig= 12
              Basis= 46,  Time=  10.50,  New conv eig= 14
              Basis= 47,  Time=  10.81,  New conv eig= 15
              Basis= 48,  Time=  11.16,  New conv eig= 16
              Basis= 49,  Time=  11.42,  New conv eig= 17
              Basis= 50,  Time=  11.73,  New conv eig= 17
              Basis= 51,  Time=  12.05,  New conv eig= 18
              Basis= 52,  Time=  12.37,  New conv eig= 18
              Basis= 53,  Time=  12.65,  New conv eig= 19
              Basis= 54,  Time=  13.01,  New conv eig= 19
              Basis= 55,  Time=  13.41,  New conv eig= 22
              Basis= 56,  Time=  13.70,  New conv eig= 24
              Basis= 57,  Time=  14.00,  New conv eig= 28
End of sweep: Basis= 57,  Time=  14.00,  New conv eig= 28
              Basis= 38,  Time=  16.12,  New conv eig=  0
              Basis= 39,  Time=  16.30,  New conv eig=  0
              Basis= 40,  Time=  16.52,  New conv eig=  0
              Basis= 41,  Time=  16.72,  New conv eig=  0
              Basis= 42,  Time=  16.94,  New conv eig=  0
End of sweep: Basis= 42,  Time=  16.94,  New conv eig=  0

ans =

    9.6481

Display the first eigenmode, and compare it to the built-in membrane plot.

pdeplot(model,'XYData',v(:,1),'ZData',v(:,1))

figure
membrane(1,20,9,9)      % the MATLAB function

Compute the sixteenth eigenvalue, and plot the sixteenth eigenmode.

l(16)                   % sixteenth eigenvalue
ans =

   92.4658

figure
pdeplot(model,'XYData',v(:,16),'ZData',v(:,16))    % sixteenth eigenmode

Compute the eigenvalues that are less than 100, and compute the corresponding eigenmodes for $- \nabla u = \lambda u$ on the geometry of the L-shaped membrane, using the legacy syntax.

Use the geometry in lshapeg. For more information about this syntax, see Create Geometry Using a Geometry Function.

g = @lshapeg;
pdegplot(g,'EdgeLabels','on')
axis equal
ylim([-1.1,1.1])

Set zero Dirichlet boundary conditions using the lshapeb function. For more information about this syntax, see Boundary Conditions by Writing Functions.

b = @lshapeb;

Set coefficients c = 1, a = 0, and d = 1. Collect eigenvalues up to 100.

c = 1;
a = 0;
d = 1;
r = [-Inf 100];

Generate a mesh and solve the eigenvalue problem.

[p,e,t] = initmesh(g,'Hmax',0.02);
[v,l] = pdeeig(b,p,e,t,c,a,d,r);
              Basis= 10,  Time=   2.89,  New conv eig=  0
              Basis= 11,  Time=   3.11,  New conv eig=  0
              Basis= 12,  Time=   3.32,  New conv eig=  0
              Basis= 13,  Time=   3.54,  New conv eig=  0
              Basis= 14,  Time=   3.73,  New conv eig=  0
              Basis= 15,  Time=   3.92,  New conv eig=  1
              Basis= 16,  Time=   4.04,  New conv eig=  1
              Basis= 17,  Time=   4.23,  New conv eig=  3
              Basis= 18,  Time=   4.42,  New conv eig=  4
              Basis= 19,  Time=   4.61,  New conv eig=  4
              Basis= 20,  Time=   4.81,  New conv eig=  4
              Basis= 21,  Time=   4.93,  New conv eig=  4
              Basis= 22,  Time=   5.12,  New conv eig=  4
              Basis= 23,  Time=   5.22,  New conv eig=  4
              Basis= 24,  Time=   5.42,  New conv eig=  5
              Basis= 25,  Time=   5.64,  New conv eig=  5
              Basis= 26,  Time=   5.83,  New conv eig=  5
              Basis= 27,  Time=   5.93,  New conv eig=  6
              Basis= 28,  Time=   6.15,  New conv eig=  7
              Basis= 29,  Time=   6.35,  New conv eig=  7
              Basis= 30,  Time=   6.52,  New conv eig=  7
              Basis= 31,  Time=   6.78,  New conv eig=  7
              Basis= 32,  Time=   7.00,  New conv eig=  8
              Basis= 33,  Time=   7.12,  New conv eig=  8
              Basis= 34,  Time=   7.44,  New conv eig=  8
              Basis= 35,  Time=   7.76,  New conv eig=  9
              Basis= 36,  Time=   7.98,  New conv eig=  9
              Basis= 37,  Time=   8.25,  New conv eig=  9
              Basis= 38,  Time=   8.55,  New conv eig=  9
              Basis= 39,  Time=   8.79,  New conv eig=  9
              Basis= 40,  Time=   9.06,  New conv eig=  9
              Basis= 41,  Time=   9.37,  New conv eig=  9
              Basis= 42,  Time=   9.73,  New conv eig= 10
              Basis= 43,  Time=   9.93,  New conv eig= 10
              Basis= 44,  Time=  10.24,  New conv eig= 12
              Basis= 45,  Time=  10.56,  New conv eig= 12
              Basis= 46,  Time=  10.85,  New conv eig= 14
              Basis= 47,  Time=  11.14,  New conv eig= 15
              Basis= 48,  Time=  11.45,  New conv eig= 16
              Basis= 49,  Time=  11.73,  New conv eig= 17
              Basis= 50,  Time=  12.02,  New conv eig= 17
              Basis= 51,  Time=  12.32,  New conv eig= 18
              Basis= 52,  Time=  12.59,  New conv eig= 18
              Basis= 53,  Time=  12.84,  New conv eig= 19
              Basis= 54,  Time=  13.17,  New conv eig= 19
              Basis= 55,  Time=  13.51,  New conv eig= 22
              Basis= 56,  Time=  13.88,  New conv eig= 24
              Basis= 57,  Time=  14.24,  New conv eig= 28
End of sweep: Basis= 57,  Time=  14.24,  New conv eig= 28
              Basis= 38,  Time=  16.42,  New conv eig=  0
              Basis= 39,  Time=  16.61,  New conv eig=  0
              Basis= 40,  Time=  16.80,  New conv eig=  0
              Basis= 41,  Time=  17.01,  New conv eig=  0
              Basis= 42,  Time=  17.24,  New conv eig=  0
End of sweep: Basis= 42,  Time=  17.24,  New conv eig=  0

Find the first eigenvalue.

l(1)
ans =

    9.6481

Import a simple 3-D geometry and find eigenvalues and eigenvectors from the associated finite element matrices.

Create a model and import the BracketWithHole.stl geometry.

model = createpde();
importGeometry(model,'BracketWithHole.stl');
figure
pdegplot(model,'FaceLabels','on')
view(30,30)
title('Bracket with Face Labels')
figure
pdegplot(model,'FaceLabels','on')
view(-134,-32)
title('Bracket with Face Labels, Rear View')

Set coefficients c = 1, a = 0, and d = 1. Collect eigenvalues that are less than 100.

c = 1;
a = 0;
d = 1;
r = [-Inf 100];

Generate a mesh for the model.

generateMesh(model);

Create the associated finite element matrices.

[Kc,~,B,~] = assempde(model,c,a,0);
[~,M,~] = assema(model,0,d,0);

Solve the eigenvalue problem.

[v,l] = pdeeig(Kc,B,M,r);
              Basis= 10,  Time=   8.79,  New conv eig=  0
              Basis= 11,  Time=   9.00,  New conv eig=  0
              Basis= 12,  Time=   9.21,  New conv eig=  0
              Basis= 13,  Time=   9.40,  New conv eig=  1
              Basis= 14,  Time=   9.60,  New conv eig=  2
              Basis= 15,  Time=   9.83,  New conv eig=  2
              Basis= 16,  Time=  10.04,  New conv eig=  3
End of sweep: Basis= 16,  Time=  10.04,  New conv eig=  3
              Basis= 13,  Time=  12.51,  New conv eig=  0
End of sweep: Basis= 13,  Time=  12.51,  New conv eig=  0

Look at the first two eigenvalues.

l([1,2])
ans =

   -0.0000
   42.8170

Plot the solution corresponding to eigenvalue 2.

pdeplot3D(model,'ColorMapData',v(:,2))

Input Arguments

collapse all

PDE model, specified as a PDEModel object.

Example: model = createpde(1)

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=λdu on Ω

or the system PDE eigenvalue problem

(cu)+au=λdu on Ω

There are a wide variety of ways of specifying c, detailed in c Coefficient for Systems. See also Specify Scalar PDE Coefficients in Character Form, Specify 2-D Scalar Coefficients in Function Form, and Specify 3-D PDE Coefficients in Function Form.

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

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

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=λdu on Ω

or the system PDE eigenvalue problem

(cu)+au=λdu on Ω

There are a wide variety of ways of specifying a, detailed in a or d Coefficient for Systems. See also Specify Scalar PDE Coefficients in Character Form, Specify 2-D Scalar Coefficients in Function Form, and Specify 3-D PDE Coefficients in Function Form.

Example: 2*eye(3)

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

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

(cu)+au=λdu on Ω

or the system PDE eigenvalue problem

(cu)+au=λdu on Ω

There are a wide variety of ways of specifying d, detailed in a or d Coefficient for Systems. See also Specify Scalar PDE Coefficients in Character Form, Specify 2-D Scalar Coefficients in Function Form, and Specify 3-D PDE Coefficients in Function Form.

Example: 2*eye(3)

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

Eigenvalue range, specified as a two-element real vector. Real parts of eigenvalues λ fall in the range r(1) ≤ λ ≤ r(2). r(1) can be -Inf. The algorithm returns all eigenvalues in this interval in the l output, up to a maximum of 99 eigenvalues.

Example: [-Inf,100]

Data Types: double

Boundary conditions, specified as a boundary matrix or boundary file. Pass a boundary file as a function handle or as a file name.

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

Data Types: double | char | function_handle

Mesh points, specified as a 2-by-Np matrix of points, where Np is the number of points in the mesh. For a description of the (p,e,t) matrices, see Mesh Data.

Typically, you use the p, e, and t data exported from the PDE app, or generated by initmesh or refinemesh.

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

Data Types: double

Mesh edges, specified as a 7-by-Ne matrix of edges, where Ne is the number of edges in the mesh. For a description of the (p,e,t) matrices, see Mesh Data.

Typically, you use the p, e, and t data exported from the PDE app, or generated by initmesh or refinemesh.

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

Data Types: double

Mesh triangles, specified as a 4-by-Nt matrix of triangles, where Nt is the number of triangles in the mesh. For a description of the (p,e,t) matrices, see Mesh Data.

Typically, you use the p, e, and t data exported from the PDE app, or generated by initmesh or refinemesh.

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

Data Types: double

Stiffness matrix, specified as a sparse matrix or as a full matrix. See Elliptic Equations. Typically, Kc is the output of assempde.

Dirichlet nullspace, returned as a sparse matrix. See Algorithms. Typically, B is the output of assempde.

Mass matrix. specified as a sparse matrix or a full matrix. See Elliptic Equations.

To obtain the input matrices for pdeeig, hyperbolic or parabolic, run both assema and assempde:

[Kc,Fc,B,ud] = assempde(model,c,a,f);
[~,M,~] = assema(model,0,d,f);

    Note:   Create the M matrix using assema with d, not a, as the argument before f.

Data Types: double
Complex Number Support: Yes

Output Arguments

collapse all

Eigenvectors, returned as a matrix. Suppose

  • Np is the number of mesh nodes

  • N is the number of equations

  • ev is the number of eigenvalues returned in l

Then v has size Np*N-by-ev. Each column of v corresponds to the eigenvectors of one eigenvalue. In each column, the first Np elements correspond to the eigenvector of equation 1 evaluated at the mesh nodes, the next Np elements correspond to equation 2, etc.

    Note:   Eigenvectors are determined only up to multiple by a scalar, including a negative scalar.

Eigenvalues, returned as a vector. The real parts of l are in the interval r. The real parts of l are monotone increasing.

Limitations

In the standard case c and d are positive in the entire region. All eigenvalues are positive, and 0 is a good choice for a lower bound of the interval. The cases where either c or d is zero are discussed next.

  • If d = 0 in a subregion, the mass matrix M becomes singular. This does not cause any trouble, provided that c > 0 everywhere. The pencil (K,M) has a set of infinite eigenvalues.

  • If c = 0 in a subregion, the stiffness matrix K becomes singular, and the pencil (K,M) has many zero eigenvalues. With an interval containing zero, pdeeig goes on for a very long time to find all the zero eigenvalues. Choose a positive lower bound away from zero but below the smallest nonzero eigenvalue.

  • If there is a region where both c = 0 and d = 0, we get a singular pencil. The whole eigenvalue problem is undetermined, and any value is equally plausible as an eigenvalue.

Some of the awkward cases are detected by pdeeig. If the shifted matrix is singular, another shift is attempted. If the matrix with the new shift is still singular a good guess is that the entire pencil (K,M) is singular.

If you try any problem not belonging to the standard case, you must use your knowledge of the original physical problem to interpret the results from the computation.

More About

collapse all

Tips

  • The equation coefficients cannot depend on the solution u or its gradient.

See Also

Introduced before R2006a

Was this topic helpful?