Documentation

This is machine translation

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

Note: This page has been translated by MathWorks. Please click here
To view all translated materals including this page, select Japan from the country navigator on the bottom of this page.

Vibration Of a Circular Membrane Using the MATLAB eigs Function

This example shows the calculation of the vibration modes of a circular membrane. The calculation of vibration modes requires the solution of the eigenvalue partial differential equation (PDE). In this example the solution of the eigenvalue problem is performed using both the PDE Toolbox™ solvepdeeig solver and the core MATLAB™ eigs eigensolver.

The main objective of this example is to show how eigs can be used with PDE Toolbox. Generally, the eigenvalues calculated by solvepdeeig and eigs are practically identical. However, sometimes, it is simply more convenient to use eigs than solvepdeeig. One example of this is when it is desired to calculate a specified number of eigenvalues in the vicinity of a user-specified target value. solvepdeeig requires that a lower and upper bound surrounding this target value be specified. eigs requires only that the target eigenvalue and the desired number of eigenvalues be specified.

Create a PDE Model

numberOfPDE = 1;
model = createpde(numberOfPDE);

Include Geometry

The geometry for a circle can easily be defined as shown below.

radius = 2;
g = decsg([1 0 0 radius]','C1',('C1')');

Create a geometry object and append it to the PDE Model.

geometryFromEdges(model,g);

Plot the geometry and display the edge labels for use in the boundary condition definition.

figure; 
pdegplot(model,'EdgeLabels','on'); 
axis equal
title 'Geometry With Edge Labels Displayed';

Define PDE Coefficients and Boundary Conditions

c = 1e2;
a = 0;
f = 0;
d = 10;
specifyCoefficients(model,'m',0,'d',d,'c',c,'a',a,'f',f);

Solution is zero at all four outer edges of the circle.

bOuter = applyBoundaryCondition(model,'dirichlet','Edge',(1:4),'u',0);

Generate Mesh

generateMesh(model,'Hmax',0.2);

Solve the Eigenvalue Problem Using eigs

Use assembleFEMatrices to calculate the global finite element mass and stiffness matrices with boundary conditions imposed using nullspace approach.

FEMatrices = assembleFEMatrices(model,'nullspace');
K = FEMatrices.Kc;
B = FEMatrices.B;
M = FEMatrices.M;
sigma = 1e2; 
numberEigenvalues = 5;
[eigenvectorsEigs,eigenvaluesEigs] = eigs(K,M,numberEigenvalues,sigma);

Reshape the diagonal eigenvaluesEigs matrix into a vector.

eigenvaluesEigs = diag(eigenvaluesEigs);

Find the largest eigenvalue and its index in the eigenvalues vector.

[maxEigenvaluesEigs, maxIndex] = max(eigenvaluesEigs);

Transform the eigenvectors with constrained equations removed to the full eigenvector including constrained equations.

eigenvectorsEigs = B*eigenvectorsEigs;

Solve the Eigenvalue Problem Using solvepdeeig

Set the range for solvepdeeig to be slightly larger than the range from eigs.

r = [min(eigenvaluesEigs)*.99 max(eigenvaluesEigs)*1.01];
result = solvepdeeig(model,r);
              Basis= 10,  Time=   0.03,  New conv eig=  0
              Basis= 14,  Time=   0.04,  New conv eig=  2
              Basis= 18,  Time=   0.04,  New conv eig=  2
              Basis= 22,  Time=   0.05,  New conv eig=  3
              Basis= 26,  Time=   0.12,  New conv eig=  5
End of sweep: Basis= 26,  Time=   0.12,  New conv eig=  5
              Basis= 15,  Time=   0.17,  New conv eig=  0
End of sweep: Basis= 15,  Time=   0.17,  New conv eig=  0
eigenvectorsPde = result.Eigenvectors;
eigenvaluesPde = result.Eigenvalues;

Compare the Solutions Computed by eigs and solvepdeeig

eigenValueDiff = sort(eigenvaluesPde) - sort(eigenvaluesEigs);
fprintf('Maximum difference in eigenvalues from solvepdeeig and eigs: %e\n', ...
  norm(eigenValueDiff,inf));
Maximum difference in eigenvalues from solvepdeeig and eigs: 1.126921e-11

As can be seen, both functions calculate the same eigenvalues. For any eigenvalue, the eigenvector can be multiplied by an arbitrary scalar. eigs and solvepdeeigs choose a different arbitrary scalar for normalizing their eigenvectors as shown in the figure below.

h = figure;
h.Position = [1 1 2 1].*h.Position;
subplot(1,2,1); 
axis equal
pdeplot(model,'XYData', eigenvectorsEigs(:,maxIndex),'Contour','on');
title(sprintf('eigs eigenvector, eigenvalue: %12.4e', eigenvaluesEigs(maxIndex)));
xlabel('x'); 
ylabel('y');
subplot(1,2,2); 
axis equal
pdeplot(model,'XYData',eigenvectorsPde(:,end),'Contour','on');
title(sprintf('solvepdeeig eigenvector, eigenvalue: %12.4e',eigenvaluesPde(end)));
xlabel('x'); 
ylabel('y');

Was this topic helpful?