Documentation |
This section describes the solution of some eigenvalue PDE problems. The problems are solved using the PDE app and command-line functions. The problems include:
On this page… |
---|
Eigenvalues and Eigenfunctions for the L-Shaped Membrane |
The problem of finding the eigenvalues and the corresponding eigenfunctions of an L-shaped membrane is of interest to all MATLAB^{®} users, since the plot of the first eigenfunction is the MathWorks^{®} logo. In fact, you can compare the eigenvalues and eigenfunctions computed by Partial Differential Equation Toolbox™ software to the ones produced by the MATLAB function membrane.
The problem is to compute all eigenmodes with eigenvalues < 100 for the eigenmode PDE problem
–Δu = λu
on the geometry of the L-shaped membrane. u = 0 on the boundary (Dirichlet condition).
With the PDE app active, check that the current mode is set to Generic Scalar. Then draw the L-shape as a polygon with corners in (0,0), (–1,0), (–1,–1), (1,–1), (1,1), and (0,1).
There is no need to define any boundary conditions for this problem since the default condition—u = 0 on the boundary—is the correct one. Therefore, you can continue to the next step: to initialize the mesh. Refine the initial mesh twice. Defining the eigenvalue PDE problem is also easy. Open the PDE Specification dialog box and select Eigenmodes. The default values for the PDE coefficients, c = 1, a = 0, d = 1, all match the problem description, so you can exit the PDE Specification dialog box by clicking the OK button.
Open the Solve Parameters dialog box by selecting Parameters from the Solve menu. The dialog box contains an edit box for entering the eigenvalue search range. The default entry is [0 100], which is just what you want.
Finally, solve the L-shaped membrane problem by clicking the = button. The solution displayed is the first eigenfunction. The value of the first (smallest) eigenvalue is also displayed. You find the number of eigenvalues on the information line at the bottom of the PDE app. You can open the Plot Selection dialog box and choose which eigenfunction to plot by selecting from a pop-up menu of the corresponding eigenvalues.
The geometry of the L-shaped membrane is described in the file lshapeg.m and the boundary conditions in the file lshapeb.m.
First, initialize the mesh and refine it twice using the command line functions at the MATLAB prompt:
[p,e,t] = initmesh('lshapeg'); [p,e,t] = refinemesh('lshapeg',p,e,t); [p,e,t] = refinemesh('lshapeg',p,e,t);
Recall the general eigenvalue PDE problem description:
$$-\nabla \cdot \left(c\nabla u\right)+au=\lambda du,$$
This means that in this case you have c = 1, a = 0, and d = 1. The syntax of pdeeig, the Partial Differential Equation Toolbox eigenvalue solver, is
[v,l] = pdeeig(b,p,e,t,c,a,d,r)
The input argument r is a two-element vector indicating the interval on the real axis where pdeeig searches for eigenvalues. Here you are looking for eigenvalues < 100, so the interval you use is [0 100].
Now you can call pdeeig and see how many eigenvalues you find:
[v,l] = pdeeig('lshapeb',p,e,t,1,0,1,[0 100]);
There are 19 eigenvalues smaller than 100. Plot the first eigenmode and compare it to the MATLAB membrane function:
pdesurf(p,t,v(:,1)) figure membrane(1,20,9,9)
membrane can produce the first 12 eigenfunctions for the L-shaped membrane. Compare also the 12th eigenmodes:
figure pdesurf(p,t,v(:,12)) figure membrane(12,20,9,9)
Looking at the following eigenmodes, you can see how the number of oscillations increases. The eigenfunctions are symmetric or antisymmetric around the diagonal from (0,0) to (1,-1), which divides the L-shaped membrane into two mirror images. In a practical computation, you could take advantage of such symmetries in the PDE problem, and solve over a region half the size. The eigenvalues of the full L-shaped membrane are the union of those of the half with Dirichlet boundary condition along the diagonal (eigenvalues 2, 4, 7, 11, 13, 16, and 17) and those with Neumann boundary condition (eigenvalues 1, 3, 5, 6, 10, 12, 14, and 15).
The eigenvalues λ_{8} and λ_{9} make up a double eigenvalue for the PDE at around 49.64. Also, the eigenvalues λ_{18} and λ_{19} make up another double eigenvalue at around 99.87. You may have gotten two different but close values. The default triangulation made by initmesh is not symmetric around the diagonal, but a symmetric grid gives a matrix with a true double eigenvalue. Each of the eigenfunctions u_{8} and u_{9} consists of three copies of eigenfunctions over the unit square, corresponding to its double second eigenvalue. You may not have obtained the zero values along a diagonal of the square—any line through the center of the square may have been computed. This shows a general fact about multiple eigenvalues for symmetric matrices; namely that any vector in the invariant subspace is equally valid as an eigenvector. The two eigenfunctions u_{8} and u_{9 }are orthogonal to each other if the dividing lines make right angles. Check your solutions for that.
Actually, the eigenvalues of the square can be computed exactly. They are
(m^{2} + n^{2})π^{2}
e.g., the double eigenvalue λ_{18} and λ_{19 }is 10π^{2}, which is pretty close to 100.
If you compute the FEM approximation with only one refinement, you would only find 16 eigenvalues, and you obtain the wrong solution to the original problem. You can of course check for this situation by computing the eigenvalues over a slightly larger range than the original problem.
You get some information from the printout in the MATLAB command window that is printed during the computation. For this problem, the algorithm computed a new set of eigenvalue approximations and tested for convergence every third step. In the output, you get the step number, the time in seconds since the start of the eigenvalue computation, and the number of converged eigenvalues with eigenvalues both inside and outside the interval counted.
Here is what MATLAB wrote:
Basis= 10, Time= 2.70, New conv eig= 0 Basis= 13, Time= 3.50, New conv eig= 0 Basis= 16, Time= 4.36, New conv eig= 0 Basis= 19, Time= 5.34, New conv eig= 1 Basis= 22, Time= 6.46, New conv eig= 2 Basis= 25, Time= 7.61, New conv eig= 3 Basis= 28, Time= 8.86, New conv eig= 3 Basis= 31, Time= 10.23, New conv eig= 5 Basis= 34, Time= 11.69, New conv eig= 5 Basis= 37, Time= 13.28, New conv eig= 7 Basis= 40, Time= 14.97, New conv eig= 8 Basis= 43, Time= 16.77, New conv eig= 9 Basis= 46, Time= 18.70, New conv eig= 11 Basis= 49, Time= 20.73, New conv eig= 11 Basis= 52, Time= 22.90, New conv eig= 13 Basis= 55, Time= 25.13, New conv eig= 14 Basis= 58, Time= 27.58, New conv eig= 14 Basis= 61, Time= 30.13, New conv eig= 15 Basis= 64, Time= 32.83, New conv eig= 16 Basis= 67, Time= 35.64, New conv eig= 18 Basis= 70, Time= 38.62, New conv eig= 22 End of sweep: Basis= 70, Time= 38.62, New conv eig= 22 Basis= 32, Time= 43.29, New conv eig= 0 Basis= 35, Time= 44.70, New conv eig= 0 Basis= 38, Time= 46.22, New conv eig= 0 Basis= 41, Time= 47.81, New conv eig= 0 Basis= 44, Time= 49.52, New conv eig= 0 Basis= 47, Time= 51.35, New conv eig= 0 Basis= 50, Time= 53.27, New conv eig= 0 Basis= 53, Time= 55.30, New conv eig= 0 End of sweep: Basis= 53, Time= 55.30, New conv eig= 0
You can see that two Arnoldi runs were made. In the first, 22 eigenvalues converged after a basis of size 70 was computed; in the second, where the vectors were orthogonalized against all the 22 converged vectors, the smallest eigenvalue stabilized at a value outside of the interval [0, 100], so the algorithm signaled convergence. Of the 22 converged eigenvalues, 19 were inside the search interval.
An extension of this problem is to compute the eigenvalues for an L-shaped membrane where the inner corner at the "knee" is rounded. The roundness is created by adding a circle so that the circle's arc is a part of the L-shaped membrane's boundary. By varying the circle's radius, the degree of roundness can be controlled. The lshapec file is an extension of an ordinary model file created using the PDE app. It contains the lines
pdepoly([-1, 1, 1, 0, 0, -1],... [-1, -1, 1, 1, 0, 0],'P1'); pdecirc(-a,a,a,'C1'); pderect([-a 0 a 0],'SQ1');
The extra circle and rectangle that are added using pdecirc and pderect to create the rounded corner are affected by the added input argument a through a couple of extra lines of MATLAB code. This is possible since Partial Differential Equation Toolbox software is a part of the open MATLAB environment.
With lshapec you can create L-shaped rounded geometries with different degrees of roundness. If you use lshapec without an input argument, a default radius of 0.5 is used. Otherwise, use lshapec(a), where a is the radius of the circle.
Experimenting using different values for the radius a shows you that the eigenvalues and the frequencies of the corresponding eigenmodes decrease as the radius increases, and the shape of the L-shaped membrane becomes more rounded. In the following figure, the first eigenmode of an L-shaped membrane with a rounded corner is plotted.
First Eigenmode for an L-Shaped Membrane with a Rounded Corner
Let us study the eigenvalues and eigenmodes of a square with an interesting set of boundary conditions. The square has corners in (-1,-1), (-1,1), (1,1), and (1,-1). The boundary conditions are as follows:
On the left boundary, the Dirichlet condition u = 0.
On the upper and lower boundary, the Neumann condition
$$\frac{\partial u}{\partial n}=0.$$
On the right boundary, the generalized Neumann condition
$$\frac{\partial u}{\partial n}-\frac{3}{4}u=0.$$
The eigenvalue PDE problem is
–Δu = λu .
We are interested in the eigenvalues smaller than 10 and the corresponding eigenmodes, so the search range is [-Inf 10]. The sign in the generalized Neumann condition is such that there are negative eigenvalues.
Using the PDE app in the generic scalar mode, draw the square using the Rectangle/square option from the Draw menu or the button with the rectangle icon. Then define the boundary conditions by clicking the ∂Ω button and then double-click the boundaries to define the boundary conditions. On the right side boundary, you have the generalized Neumann conditions, and you enter them as constants: g = 0 and g = –3/4.
Initialize the mesh and refine it once by clicking the Δ and refine buttons or by selecting the corresponding options from the Mesh menu.
Also, define the eigenvalue PDE problem by opening the PDE Specification dialog box and selecting the Eigenmodes option. The general eigenvalue PDE is described by
$$-\nabla \cdot \left(c\nabla u\right)+au=\lambda du,$$
so for this problem you use the default values c = 1, a = 0, and d = 1. Also, in the Solve Parameters dialog box, enter the eigenvalue range as the MATLAB vector [-Inf 10].
Finally, click the = button to compute the solution. By default, the first eigenfunction is plotted. You can plot the other eigenfunctions by selecting the corresponding eigenvalue from a pop-up menu in the Plot Selection dialog box. The pop-up menu contains all the eigenvalues found in the specified range. You can also export the eigenfunctions and eigenvalues to the MATLAB main workspace by using the Export Solution option from the Solve menu.
This example shows how to compute the eigenvalues and eigenmodes of a square domain using command-line functions. The geometry description file and boundary condition file for this problem are called squareg.m and squareb2.m, respectively. Create and refine the mesh for the problem:
[p,e,t]=initmesh('squareg'); [p,e,t]=refinemesh('squareg',p,e,t);
The eigenvalue PDE coefficients for this problem are c = 1, a = 0, and d = 1. You can enter the eigenvalue range r as the vector [-Inf 10]. pdeeig returns two output arguments, the eigenvalues as an array l and a matrix v of corresponding eigenfunctions:
[v,l]=pdeeig('squareb2',p,e,t,1,0,1,[-Inf 10]);
Basis= 10, Time= 0.15, New conv eig= 0 Basis= 17, Time= 0.16, New conv eig= 2 Basis= 24, Time= 0.17, New conv eig= 7 End of sweep: Basis= 24, Time= 0.18, New conv eig= 7 Basis= 17, Time= 0.27, New conv eig= 0 End of sweep: Basis= 17, Time= 0.27, New conv eig= 0
To plot the fourth eigenfunction as a surface plot, enter
pdesurf(p,t,v(:,4))
This problem is separable, meaning
The functions f and g are eigenfunctions in the x and y directions, respectively. In the x direction, the first eigenmode is a slowly increasing exponential function. The higher modes include sinusoids. In the y direction, the first eigenmode is a straight line (constant), the second is half a cosine, the third is a full cosine, the fourth is one and a half full cosines, etc. These eigenmodes in the y direction are associated with the eigenvalues
There are five eigenvalues smaller than 10 for this problem, and the first one is even negative (-0.4145). It is possible to trace the preceding eigenvalues in the eigenvalues of the solution. Looking at a plot of the first eigenmode, you can see that it is made up of the first eigenmodes in the x and y directions. The second eigenmode is made up of the first eigenmode in the x direction and the second eigenmode in the y direction.
Look at the difference between the first and the second eigenvalue compared to :
l(2)-l(1)
ans = 2.4745
pi^2/4
ans = 2.4674
Likewise, the fifth eigenmode is made up of the first eigenmode in the x direction and the third eigenmode in the y direction. As expected, l(5)-l(1) is approximately equal to :
l(5) - l(1) - pi^2
ans = 0.0685
You can explore higher modes by increasing the search range to include eigenvalues greater than 10.