Documentation

## Vibration of Circular Membrane

This example shows how to calculate the vibration modes of a circular membrane by using the MATLAB `eigs` function.

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's `solvepdeeig` solver and the core MATLAB's `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);```

Create a geometry object and include it in the PDE model. The geometry is a circle.

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

Plot the geometry and display the edge labels.

```figure; pdegplot(model,'EdgeLabels','on'); axis equal title 'Geometry With Edge Labels Displayed';``` Specify the coefficients.

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

Specify the boundary conditions: the solution is zero at all four outer edges of the circle.

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

Generate a mesh.

`generateMesh(model,'Hmax',0.2);`

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;```

Solve the eigenvalue problem by using the `eigs` function.

```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 by using the `solvepdeeig` function. 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.19, New conv eig= 0 Basis= 13, Time= 0.27, New conv eig= 2 Basis= 16, Time= 0.33, New conv eig= 2 Basis= 19, Time= 0.35, New conv eig= 2 Basis= 22, Time= 0.39, New conv eig= 3 Basis= 25, Time= 0.45, New conv eig= 3 Basis= 28, Time= 0.46, New conv eig= 5 End of sweep: Basis= 28, Time= 0.46, New conv eig= 5 Basis= 15, Time= 0.48, New conv eig= 0 End of sweep: Basis= 15, Time= 0.48, 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: 9.947598e-14 ```

As you can see, both functions calculate the same eigenvalues. For any eigenvalue, the eigenvector can be multiplied by an arbitrary scalar. The `eigs` and `solvepdeeig` functions 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');``` ##### Support Get trial now