LAB 12: NORMAL MODES OF A MEMBRANE AND BOUNDARY-VALUE ODE PROBLEMS

Mathematics:

Normal modes of mechanical vibrations of a circular membrane can be found by solving the boundary-value problems for a second-order ordinary differential equation. General vibrations of a circular membrane are decomposed into a superposition (linear combination) of individual normal modes. Each mode has its own shape and oscillates with its own frequency. This project exploits radially symmetrical oscillations of a circular membrane with a hole at the center.

Denote r as the radius from the membrane's center in polar coordinates on the plane: r =. The membrane is located between: a  r  b and has a hole at r < a. The shape of the radially symmetric normal mode of the membrane is found from the second-order differential equation:

+  + y = 0                                                                (1)

with the boundary conditions at the boundaries of the membrane:

y(a) = 0,    y(b) = 0                                                                                (2)

The parameter  is the frequency of oscillations of the normal mode. If u = u(r,t) is the vertical displacement of the membrane as function of time t, then: u(r,t) = y(r) cos(t). The problem (1)-(2) is a boundary-value problem for y = y(r), since the two conditions that define a solution y = y(r) are given at the two different ends of the membrane: r = a and r = b. It is clear from (1)-(2) that y(r) = 0 is always a solution (the trivial zero solution). However, boundary-value problems may not have a non-trivial solution such that y(r)  0. A non-trivial non-zero solution for y = y(r) exists only for a specific value of the parameter . The parameter  is called the eigenvalue and the function y(r) is called the eigenfunction of the boundary-value problem (1)-(2).

There are several possible non-zero solutions of the boundary-value problem (1)-(2). Let us enumerate these solutions by a sequence{  } such that  <  <   < … with the corresponding eigenfunctions { yk(r) }. The eigenfunctions are classified by the number of zeros (nodes) in between a < r < b. The first (smallest) eigenvalue  corresponds to eigenfunction y1(r) with no nodes (zeros). The second eigenvalue  corresponds to the eigenfunction y2(r) with one node and so on. The eigenvalue  corresponds to eigenfunction yk(r) with exactly (k-1) nodes between a <  r < b .

If the circular membrane has no hole, a = 0 and the eigenfunctions of the boundary-value problem (1)-(2) coincide with the Bessel function of zero order: yk(r) = J0(r) while the eigenvalues are given by zeros of the Bessel function at r = b: J0(b) = 0. The first five eigenvalues  are shown in the table:

 eigenvalues k = 1 k = 2 k = 3 k = 4 k = 5 b 2.40482555770 5.2007811029 8.65372791291 11.79153443901 14.93091770849

Objectives:

·         Understand how to solve the boundary-value problem for a second-order differential equation with the central difference method

·         Find the first five eigenvalues and eigenfunctions of the circular membrane

·         Exploit MATLAB solvers for boundary-value ODE problems

#### MATLAB script for finding eigenvalues of the circular membrane with the central difference method

In numerical methods of solving the boundary-value problem for a second-order differential equation, the solution is to be found on a discrete grid: r1,r2,…,rn,rn+1, where r1 = a and rn+1 = b. If the grid is equally spaced with a constant step size h, then h = (b-a) / n. Derivatives of a solution y = y(r) at any point of the grid r = rk can be approximated with the central differences:

(rk) =      (rk)  =

This is the central difference method for solving the boundary-value problem for a second-order differential equation. Replacing the derivatives by the central differences at each interior point of the grid: r = rk for k = 2,…,n; the second-order differential equation becomes a linear system of equations:

yk+1 ( 1 + ) + yk ( -2 + h2 ) + yk-1 ( 1 - ) = 0

With the boundary values y1 = yn+1 = 0, the linear system is converted into the matrix form:

M y =  y

where  = h2 and the vector y and the coefficient matrix M is for the case n = 5:

y =       M =

This is a linear eigenvalue problem for the matrix M with the eigenvalues  = h2. First five eigenvalues of the matrix M are computed and summarized in the table:

 eigenvalue k = 1 k = 2 k = 3 k = 4 k = 5 0.7630790233 1.5556700046 2.3412414511 3.1214853270 3.8961027211

Steps in writing the MATLAB script:

1. Define a = 1, b = 5. Set up a partition for the interval a < r < b with constant step size h = 0.1. Find n.
2. Set up the matrix M for the boundary-value problem.
3. Find eigenvalues and eigenvectors of the matrix M with MATLAB function "eig".
4. Find frequenciesfrom the eigenvalues of the matrix M. Output the first five smallest eigenvalues.
5. Plot the eigenvectors corresponding to the first five eigenvalues . Check that the number of zeros of the k-th eigenvector between a < r < b equals exactly (k-1).

Exploiting the MATLAB script:

1. Repeat computations with h = 0.05 and 0.025. Compare the values of the first five eigenvalues.
2. Use Richardson extrapolation twice to improve the numerical approximation for the first eigenvalue.

If the central difference approximation gives the order O(h2) truncation error to the eigenvalue = D1(h), then the Richardson extrapolations predict the order O(h4) and O(h6) truncation errors to the approximations:

D2(h) =      D3(h) =

#### MATLAB script for computing eigenfunctions with MATLAB ODE solvers

The eigenfunction of the circular membrane y = y(r) solves the following initial-value problem:

= -  - y

with the initial conditions at the inner boundary of the membrane:

y(a) = 0,    y'(a) = 1

The slope of y = y(r) at r = a is normalized to one (any solution y = y(r) can be multiplied by a non-zero constant). The solution of the initial-value problem y = y(r) satisfies also the boundary condition at the outer boundary of the membrane y(b) = 0 if   = is the eigenvalue of the circular membrane. If   is close but not equal to, the eigenfunction y = y(r) is close but not equal to zero at r = b.

The second-order differential equation for y = y(r) can be rewritten as a system of two first-order differential equations:

y' = u                                                         (1)

u' = -  u –  y                                     (2)

with the initial conditions:

y(a) = 0,     u(a) = 1                              (3)

The initial-value problem can be solved if the value  is known explicitly or it is approximated by results of the central difference method.

Steps in writing the MATLAB script:

1. Define the system of differential equations (1)-(2) in a separate function M-file "system" with three input arguments: r as scalar, Y = [y,u] as vector, and omega as a parameter of the system, and with one output argument Ydot as vector of right-hand-sides of the system.
2. Define a vector r for a partition of the interval [a, b] with step size h = 0.1.
3. Define a vector for initial values: Y0 = [0,1].
4. Define parameter  according to the central difference approximations for the first five eigenvalues. For each eigenvalue, call MATLAB ODE solver "ode23" to solve the system of differential equations (1)-(2).
5. Retrive the output of the ODE solver "ode23" and plot y as a function of r for first five eigenvectors.

Exploiting the MATLAB script:

1. Compare the eigenvectors obtained by the central difference method and by the ODE solver "ode23".
2. Estimate the error of numerical approximation for  by comparing the value y(b) and 0 in the solution of the system of differential equations.

#### Compute the first five eigenvalues (frequencies) for non-symmetric normal modes of the circular membrane:

+  - y + y = 0,      n = 1,2

Output the first five eigenvalues for n = 1 and n = 2 separately.