(*:Name: Chladni` *) (*:Author: Wilberd van der Kallen, 2006 *) (*:Mathematica Version: 5.2 *) (*:Package Version: 1.1 *) (*:Summary: *) (* We introduce functions that model some Chladni plate experiments. *) (*:Sources: *) (* Wheatstone, On the Figures Obtained by Strewing Sand on Vibrating *) (* Surfaces, Commonly Called Acoustic Figures, Philosophical Transactions *) (* of the Royal Society of London, vol 123, 1833. *) (* Lord Rayleigh, Theory of Sound, Dover, New York 1945. *) (* Mary Waller, Chladni Figures, A Study in Symmetry, Bell, London 1961. *) (* Rossing en Fletcher, Principles of Vibration and Sound, Springer, *) (* New York 1995. *) (* :Discussion: *) (* Consider a snapshot of an oscillating thin bending bar with free ends, *) (* in an eigen mode. *) (* The equations are *) (* u''''= freq^2 u *) (* with boundary conditions *) (* u'''(1) = u'''(-1) = u''(1) = u''(-1) = 0. *) (* Here freq is proportional to the frequency of the oscillation. *) (* The proportionality constant depends on things like the *) (* stiffness of the bar. *) (* *) (* Next consider a snapshot for an oscillating thin *) (* round plate of radius one with free boundary, in an eigen mode. *) (* We simplify the boundary conditions. *) (* If Delta denotes the Laplacian, then the equations are *) (* Delta Delta u = freq^2 u *) (* with our simplified boundary conditions saying that both *) (* Delta u *) (* and its gradient vanish on the boundary. *) (* One gets a model that, given its simplicity, agrees reasonably well *) (* with reality, both in shape of Chladni pattern and in frequency. *) (* *) (* Next consider a snapshot of an oscillating thin bending square plate *) (* with free boundary. It is too hard to solve the equations *) (* symbolically. By way of approximation one could work with the pretend *) (* solution f(x)g(y) + (or -) f(y)g(x), where f and g are solutions *) (* for the bending bar. This pretend solution does not solve any *) (* equation. If f and g satisfy the free boundary conditions for the *) (* bending bar, then the product f(x)g(y) is no eigenfunction of *) (* Delta Delta. But if one removes the exponentials from f and g, *) (* keeping only the sines and cosines, then one does get an eigenfunction *) (* of Delta Delta, with a frequency that is the sum of the frequencies *) (* associated to f and g. Thus one gets a reasonable choice of frequency *) (* to associate with the pretend solution. *) (* One gets a very rough model that shows some similarity with reality, *) (* both in shape of Chladni pattern and in frequency. This goes back to *) (* Wheatstone. *) BeginPackage["Chladni`"] bar::usage = "bar[m] gives a snapshot of an oscillating thin \ bending bar with free ends at -1 and 1." freqbar::usage = "freqbar[m] is proportional with the frequency associated with \ bar[m]." freqsquare::usage = "freqsquare[m,n] = freqbar[m] + freqbar[n]" disc::usage = "disc[m,n][r] gives the radial part of \ a snapshot of an oscillating thin disc of \ radius one with simplified free boundary. The angular part is Cos[m phi]." freqdisc::usage = "freqdisc[m,n] is proportional with the \ frequency associated with disc[m,n]." (*:Examples: In[1]:= < "disc(" <> ToString[3] <> ") freq=" <> ToString[freqbar[3]]] Out[2]= -Graphics- In[3]:= ContourPlot[ bar[7][x]bar[4][y] + bar[4][x]bar[7][y], {x, -1, 1}, {y, -1, 1}, Contours -> {-0.01, 0.01}, PlotPoints -> 100, PlotLabel -> "(" <> ToString[7] <> "," <> ToString[4] <> ")plus(" <> ToString[4] <> "," <> ToString[7] <> ") freq=" <> ToString[freqsquare[7, 4]]] Out[3]= -ContourGraphics- In[4]:= ContourPlot[(disc[3, 2][r] /. r -> Sqrt[u^2 + v^2] // Evaluate )Cos[ 3 Arg[u + I v]]If[u^2 + v^2 < 1, 1, 0], {u, -1, 1}, {v, -1, 1}, Contours -> {-.01, .01}, PlotPoints -> 100, PlotLabel -> "(" <> ToString[3] <> "," <> ToString[2] <> ") freq=" <> ToString[freqdisc[3, 2]]] Out[4]= -ContourGraphics- *) Begin["Chladni`Private`"] {a,m,b,n,x,c,nn,d,k,ka,y0,la,r}; Clear[a,m,b,n,x,c,nn,d,k,ka,y0,la,r]; bar[0] = Function[x, 1] bar[1] = Function[x, x] bar[m_Integer] := (bar[m] = If[OddQ[m], Function[x, Evaluate[First[u[x, 0, b, 0, d, k] /. threecond] /. d -> 1 /. k -> onevk[(m - 1)/2]]], Function[x, Evaluate[First[u[x, a, 0, c, 0, k] /. threecondeven] /. c -> 1 /. k -> evk[m/2]]]]) /; m > 1 u[x_, a_, b_, c_, d_, k_] = c*Cos[k*x] + a*Cosh[k*x] + d*Sin[k*x] + b*Sinh[k*x] threecond = {{a -> (d*Sin[k]*(-Sin[k] + Cos[k]*Tanh[k]))/ (Cosh[k]*Sin[k] + Cos[k]*Sinh[k]), b -> d*Cos[k]*Sech[k], c -> (d*Sinh[k]*(Sin[k] - Cos[k]*Tanh[k]))/(Cosh[k]*Sin[k] + Cos[k]*Sinh[k])}} onevk[0] = 0 onevk[m_Integer] := (onevk[m] = k /. FindRoot[trans[k] == 0, {k, onevk[m - 1] + Pi}]) /; m >= 1 trans[k_] = -2*Sin[k] + 2*Cos[k]*Tanh[k] threecondeven = {{a -> -(c*Csch[k]*Sin[k]), b -> -((c*Cos[k]*(Cos[k] + Coth[k]*Sin[k]))/(-(Cosh[k]*Sin[k]) + Cos[k]*Sinh[k])), d -> (c*Cosh[k]*(Cos[k] + Coth[k]*Sin[k]))/ (Cosh[k]*Sin[k] - Cos[k]*Sinh[k])}} evk[0] = -1 evk[m_Integer] := (evk[m] = k /. FindRoot[transeven[k] == 0, {k, evk[m - 1] + Pi}]) /; m >= 1 transeven[k_] = -2*(Cos[k] + Coth[k]*Sin[k]) freqsquare[m_, n_] = freqbar[m] + freqbar[n] freqbar[0] = 0 freqbar[1] = 0 freqbar[m_Integer] := If[OddQ[m], onevk[(m - 1)/2]^2, evk[m/2]^2] /; m > 1 freqdisc[m_Integer, 1] := (freqdisc[m, 1] = (ka /. FindRoot[err[ka, m], {ka, (38*m/100 + 1)*Pi}])^2) /; Inequality[0, LessEqual, m, Less, 11] freqdisc[m_Integer, 1] := (freqdisc[m, 1] = (ka /. FindRoot[err[ka, m], {ka, 1 + Sqrt[freqdisc[m - 1, 1]]}])^ 2) /; m > 10 freqdisc[m_Integer, 2] := (freqdisc[m, 2] = (ka /. FindRoot[err[ka, m], {ka, 1 + Sqrt[freqdisc[m - 1, 2]]}])^ 2) /; m > 5 freqdisc[m_Integer, n_Integer] := (freqdisc[m, n] = (ka /. FindRoot[err[ka, m], {ka, Sqrt[freqdisc[m, n - 1]] + Pi}])^ 2) /; m >= 0 && n > 0 err[ka_, 0] = ka*BesselJ[1, ka] - (ka^2*(BesselJ[0, ka] - BesselJ[2, ka]))/ 2 - (ka*BesselI[1, ka]*(ka*BesselJ[0, ka] + 2*BesselJ[1, ka] - ka*BesselJ[2, ka]))/(ka*BesselI[0, ka] + 2*BesselI[1, ka] + ka*BesselI[2, ka]) + (ka^2*(BesselI[0, ka] + BesselI[2, ka])* (ka*BesselJ[0, ka] + 2*BesselJ[1, ka] - ka*BesselJ[2, ka]))/ (2*(ka*BesselI[0, ka] + 2*BesselI[1, ka] + ka*BesselI[2, ka])) + (ka^2*(ka*BesselI[1, ka] + (ka*(BesselI[1, ka] + BesselI[3, ka]))/2)* (ka*BesselJ[0, ka] + 2*BesselJ[1, ka] - ka*BesselJ[2, ka]))/ (2*(ka*BesselI[0, ka] + 2*BesselI[1, ka] + ka*BesselI[2, ka])) - (ka^2*(-(ka*BesselJ[1, ka]) - (ka*(BesselJ[1, ka] - BesselJ[3, ka]))/2))/ 2 err[ka_, nn_] = -(ka*(BesselJ[-1 + nn, ka] - BesselJ[1 + nn, ka]))/2 + (ka*((ka*(BesselJ[-2 + nn, ka] - BesselJ[nn, ka]))/2 - (ka*(BesselJ[nn, ka] - BesselJ[2 + nn, ka]))/2))/2 - (ka*(BesselI[-1 + nn, ka] + BesselI[1 + nn, ka])* (ka^2*BesselJ[-2 + nn, ka] + 2*ka*BesselJ[-1 + nn, ka] - 2*ka^2*BesselJ[nn, ka] - 4*nn^2*BesselJ[nn, ka] - 2*ka*BesselJ[1 + nn, ka] + ka^2*BesselJ[2 + nn, ka]))/ (2*(-(ka^2*BesselI[-2 + nn, ka]) - 2*ka*BesselI[-1 + nn, ka] - 2*ka^2*BesselI[nn, ka] + 4*nn^2*BesselI[nn, ka] - 2*ka*BesselI[1 + nn, ka] - ka^2*BesselI[2 + nn, ka])) + (ka*((ka*(BesselI[-2 + nn, ka] + BesselI[nn, ka]))/2 + (ka*(BesselI[nn, ka] + BesselI[2 + nn, ka]))/2)* (ka^2*BesselJ[-2 + nn, ka] + 2*ka*BesselJ[-1 + nn, ka] - 2*ka^2*BesselJ[nn, ka] - 4*nn^2*BesselJ[nn, ka] - 2*ka*BesselJ[1 + nn, ka] + ka^2*BesselJ[2 + nn, ka]))/ (2*(-(ka^2*BesselI[-2 + nn, ka]) - 2*ka*BesselI[-1 + nn, ka] - 2*ka^2*BesselI[nn, ka] + 4*nn^2*BesselI[nn, ka] - 2*ka*BesselI[1 + nn, ka] - ka^2*BesselI[2 + nn, ka])) + (ka*((ka*((ka*(BesselI[-3 + nn, ka] + BesselI[-1 + nn, ka]))/2 + (ka*(BesselI[-1 + nn, ka] + BesselI[1 + nn, ka]))/2))/2 + (ka*((ka*(BesselI[-1 + nn, ka] + BesselI[1 + nn, ka]))/2 + (ka*(BesselI[1 + nn, ka] + BesselI[3 + nn, ka]))/2))/2)* (ka^2*BesselJ[-2 + nn, ka] + 2*ka*BesselJ[-1 + nn, ka] - 2*ka^2*BesselJ[nn, ka] - 4*nn^2*BesselJ[nn, ka] - 2*ka*BesselJ[1 + nn, ka] + ka^2*BesselJ[2 + nn, ka]))/ (2*(-(ka^2*BesselI[-2 + nn, ka]) - 2*ka*BesselI[-1 + nn, ka] - 2*ka^2*BesselI[nn, ka] + 4*nn^2*BesselI[nn, ka] - 2*ka*BesselI[1 + nn, ka] - ka^2*BesselI[2 + nn, ka])) + 2*nn^2*(BesselJ[nn, ka] + (BesselI[nn, ka]*(ka^2*BesselJ[-2 + nn, ka] + 2*ka*BesselJ[-1 + nn, ka] - 2*ka^2*BesselJ[nn, ka] - 4*nn^2*BesselJ[nn, ka] - 2*ka*BesselJ[1 + nn, ka] + ka^2*BesselJ[2 + nn, ka]))/(-(ka^2*BesselI[-2 + nn, ka]) - 2*ka*BesselI[-1 + nn, ka] - 2*ka^2*BesselI[nn, ka] + 4*nn^2*BesselI[nn, ka] - 2*ka*BesselI[1 + nn, ka] - ka^2*BesselI[2 + nn, ka])) - nn^2*((ka*(BesselJ[-1 + nn, ka] - BesselJ[1 + nn, ka]))/2 + (ka*(BesselI[-1 + nn, ka] + BesselI[1 + nn, ka])* (ka^2*BesselJ[-2 + nn, ka] + 2*ka*BesselJ[-1 + nn, ka] - 2*ka^2*BesselJ[nn, ka] - 4*nn^2*BesselJ[nn, ka] - 2*ka*BesselJ[1 + nn, ka] + ka^2*BesselJ[2 + nn, ka]))/ (2*(-(ka^2*BesselI[-2 + nn, ka]) - 2*ka*BesselI[-1 + nn, ka] - 2*ka^2*BesselI[nn, ka] + 4*nn^2*BesselI[nn, ka] - 2*ka*BesselI[1 + nn, ka] - ka^2*BesselI[2 + nn, ka]))) + (ka*((ka*((ka*(BesselJ[-3 + nn, ka] - BesselJ[-1 + nn, ka]))/2 - (ka*(BesselJ[-1 + nn, ka] - BesselJ[1 + nn, ka]))/2))/2 - (ka*((ka*(BesselJ[-1 + nn, ka] - BesselJ[1 + nn, ka]))/2 - (ka*(BesselJ[1 + nn, ka] - BesselJ[3 + nn, ka]))/2))/2))/2 disc[d_Integer, m_Integer] := (disc[d, m] = Function[x,(la = freqdisc[d, m]^2; If[d == 0, y0[x], y[x] /. nn -> d])//Evaluate]) /; d >= 0 && m > 0 y0[r_] = BesselJ[0, la^(1/4)*r] + (BesselI[0, la^(1/4)*r]* (la^(1/4)*BesselJ[0, la^(1/4)] + 2*BesselJ[1, la^(1/4)] - la^(1/4)*BesselJ[2, la^(1/4)]))/(la^(1/4)*BesselI[0, la^(1/4)] + 2*BesselI[1, la^(1/4)] + la^(1/4)*BesselI[2, la^(1/4)]) y[r_, 0, a_, b_, la_] = a*BesselI[0, la^(1/4)*r] + b*BesselJ[0, la^(1/4)*r] y[r_] = BesselJ[nn, la^(1/4)*r] + (BesselI[nn, la^(1/4)*r]* (Sqrt[la]*BesselJ[-2 + nn, la^(1/4)] + 2*la^(1/4)* BesselJ[-1 + nn, la^(1/4)] - 2*Sqrt[la]*BesselJ[nn, la^(1/4)] - 4*nn^2*BesselJ[nn, la^(1/4)] - 2*la^(1/4)*BesselJ[1 + nn, la^(1/4)] + Sqrt[la]*BesselJ[2 + nn, la^(1/4)]))/ (-(Sqrt[la]*BesselI[-2 + nn, la^(1/4)]) - 2*la^(1/4)*BesselI[-1 + nn, la^(1/4)] - 2*Sqrt[la]*BesselI[nn, la^(1/4)] + 4*nn^2*BesselI[nn, la^(1/4)] - 2*la^(1/4)*BesselI[1 + nn, la^(1/4)] - Sqrt[la]*BesselI[2 + nn, la^(1/4)]) y[r_, n_, a_, b_, c_, d_, la_] = a*BesselI[-n, la^(1/4)*r] + c*BesselI[n, la^(1/4)*r] + b*BesselJ[-n, la^(1/4)*r] + d*BesselJ[n, la^(1/4)*r] End[] EndPackage[] Null