Documentation Center

  • Trial Software
  • Product Updates

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™ pdeeig 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 pdeeig and eigs are practically identical. However, sometimes, it is simply more convenient to use eigs than pdeeig. 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. pdeeig 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.

Geometry And Mesh

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

radius = 2;
g = decsg([1 0 0 radius]', 'C1', ('C1')');
[p,e,t] = initmesh(g, 'hmax', .2);

Define the PDE Coefficients and Boundary Conditions

c = 1e2;
a = 0;
f = 0;
d = 10;
% Define boundary conditions using a boundary file function.
b = @boundaryFileZeroDirichlet;
% This boundary file function sets the solution to zero at r=radius.
type boundaryFileZeroDirichlet.m
function [ q, g, h, r ] = boundaryFileZeroDirichlet( p, e, u, time )
%BOUNDARYFILEZERODIRICHLET Solution is zero on all edges
%   Define a Dirichlet boundary condition making the solution equal zero
%   on all boundary edges.
N = 1; % Only a single scalar PDE
ne = size(e,2); % number of boundary edges
q = zeros(N^2, ne); % Neumann coefficient q is zero on all edges
g = zeros(N, ne); % Neumann coefficient g is zero on all edges
h = ones(N^2, 2*ne); % Dirichlet h coefficient is one at both ends of all edges
r = zeros(N,2*ne); % Dirichlet r coefficient is zero at both ends of all edges
end

Solve the eigenvalue problem using eigs

Use assempde and assema to calculate the global finite element mass and stiffness matrices.

[K,~,B] = assempde(b,p,e,t,c,a,f);
[~,M] = assema(p,t,c,d,f);
M = B'*M*B; % apply the constraints to the mass matrix from |assema|
sigma = 1e2; numberEigenvalues = 5;
[eigenvectorsEigs,eigenvaluesEigs] = eigs(K,M,numberEigenvalues,sigma);
% eigs orders the eigenvalues (and their eigenvectors) from highest to
% lowest. Reorder these from lowest to highest to be consistent with |pdeeig|.
eigenvaluesEigs = flipud(diag(eigenvaluesEigs));
% Reorder the eigenvectors. Also transform the eigenvectors with constrained
% equations removed to the full eigenvector including constrained equations.
eigenvectorsEigs = B*fliplr(eigenvectorsEigs);

Solve the eigenvalue problem using pdeeig

Define the eigenvalue range for pdeeig from the output eigenvalues from eigs so that it computes the same ones.

r = [eigenvaluesEigs(1)*.99 eigenvaluesEigs(end)*1.01];
[eigenvectorsPde,eigenvaluesPde] = pdeeig(b,p,e,t,c,a,d,r);
              Basis= 10,  Time=   0.03,  New conv eig=  1
              Basis= 19,  Time=   0.03,  New conv eig=  3
              Basis= 28,  Time=   0.04,  New conv eig=  8
              Basis= 37,  Time=   0.05,  New conv eig= 12
End of sweep: Basis= 37,  Time=   0.05,  New conv eig= 12
              Basis= 22,  Time=   0.06,  New conv eig=  0
              Basis= 31,  Time=   0.07,  New conv eig=  0
End of sweep: Basis= 31,  Time=   0.07,  New conv eig=  0

Compare the solutions computed by eigs and pdeeig

eigenValueDiff = eigenvaluesPde - eigenvaluesEigs;
fprintf('Maximum difference in eigenvalues from pdeeig and eigs: %e\n', ...
  norm(eigenValueDiff,inf));
%
% As can be seen, both functions calculate the same eigenvalues. For any
% eigenvalue, the eigenvector can be multiplied by an arbitrary scalar.
% eigs and pdeeigs choose a different arbitrary scalar for normalizing
% their eigenvectors as shown in the figure below.
%
h = figure; pos = get(h,'position'); set(h,'position',[1 1 2 1].*pos);
subplot(1,2,1); axis equal;
pdeplot(p,e,t,'xydata', eigenvectorsEigs(:,end), 'contour', 'on');
title(sprintf('eigs eigenvector, eigenvalue: %12.4e', eigenvaluesEigs(end)));
xlabel('x'); ylabel('y');
subplot(1,2,2); axis equal;
pdeplot(p,e,t,'xydata', eigenvectorsPde(:,end), 'contour', 'on');
title(sprintf('pdeeig eigenvector, eigenvalue: %12.4e', eigenvaluesPde(end)));
xlabel('x'); ylabel('y');
Maximum difference in eigenvalues from pdeeig and eigs: 1.136868e-13

Was this topic helpful?