function Drum
The axisymmetric harmonic vibrations of a circular drum can be described by the ODE
$$ u''(r) + r^{-1} u'(r) = -\omega^2 u(r),~~ u'(0)=1, u(1)=0, $$
where $r$ is the radial coordinate and $\omega$ is the frequency of vibration. Only discrete positive values of $\omega$ are possible, corresponding to the eigenvalues of the differential equation.
We multiply the ODE through by $r$ to avoid a potential division by zero. This creates a generalized problem in the form $Au = \lambda Bu$.
r = chebfun('r',[0,1]); A = chebop(0,1); A.op = @(r,u) r*diff(u,2) + diff(u); A.lbc = 'neumann'; A.rbc = 'dirichlet'; B = chebop(0,1); B.op = @(r,u) r*u;
Then we find the eigenvalues with eigs
. It turns out that the $\omega$ values are also zeros of the Bessel function $J_0$, which gives a way to valudate the results.
[V,D] = eigs(A,B); [omega,ii] = sort(sqrt(-diag(D))); omega V = [V{:,ii'}]; err = omega - sort(roots( besselj(0,chebfun('r',[0 20])) ))
omega = 2.404825557946273 5.520078110504802 8.653727912948732 11.791534439018191 14.930917708488634 18.071063967910305 err = 1.0e-09 * 0.250494291975656 0.218500773030428 0.037726266555183 0.003904432333002 0.000854427639752 -0.000625277607469
We also get the eigenfunctions, which gives a way to visualize deflections of the drums for pure frequencies.
V = V*diag(sign(V(0,:))); % ensure V(0,k) > 0 [rr,tt] = meshgrid(linspace(0,1,40),linspace(0,2*pi,60)); for k = 1:4, subplot(2,2,k), mesh(rr.*cos(tt),rr.*sin(tt),repmat(V(rr(1,:),k),60,1)) zlim([-1 3]),caxis([-3 3]), view(-33,20), axis off end
If the drum instead has a variable density given by $\rho(r)$, the right- hand side of the original ODE becomes $-\omega^2\rho u$. In general, the connection to Bessel functions is broken, but we will not miss a beat using eigs
.
Constant density gives $\omega_2/\omega_1 = 2.2954$. Let's design a density so that $\omega_2/\omega_1 = 2$, a perfect octave. We will search among density functions of the form
$$ \rho(r) = 1 - a\sin(\pi r). $$
Here is a function that returns the ratio for any $a$.
function ratio = evratio(a) rho = 1 - a*sin(pi*r); B.op = @(r,u) r*rho*u; [V,D] = eigs(A,B,2,0); [omega,ii] = sort(sqrt(-diag(D))); V = [V{:,ii'}]; ratio = omega(2)/omega(1); end
Now, we create a chebfun to hit the target.
ratfun = chebfun(@evratio,[0.5,1]); astar = find(ratfun==2) clf, plot(ratfun), title('Eigenvalue ratio'), xlabel('a') set(gca,'xtick',[0.5,astar,1],'ytick',[2],'xgrid','on','ygrid','on')
astar = 0.812158808552563
We compute the ratio at astar
to verify the answer, and plot the eigenfunctions.
residual = evratio(astar) - 2 subplot(1,2,1), surfl(rr.*cos(tt),rr.*sin(tt),repmat(-V(rr(1,:),1),60,1)) shading interp, lighting phong, title('First mode') subplot(1,2,2), surfl(rr.*cos(tt),rr.*sin(tt),repmat(-V(rr(1,:),2),60,1)) shading interp, lighting phong, title('Second mode') colormap copper
residual = -3.648192858918264e-12
end
Couldn't create JOGL canvas--using painters Couldn't create JOGL canvas--using painters Couldn't create JOGL canvas--using painters Couldn't create JOGL canvas--using painters Couldn't create JOGL canvas--using painters Couldn't create JOGL canvas--using painters Couldn't create JOGL canvas--using painters Couldn't create JOGL canvas--using painters Couldn't create JOGL canvas--using painters Couldn't create JOGL canvas--using painters Couldn't create JOGL canvas--using painters Couldn't create JOGL canvas--using painters Couldn't create JOGL canvas--using painters