Staaf

Eerst de trillende dunne balk van lengte twee. 
Twee vrije uiteinden. Dat lijkt vreemd, maar is de geschikte voorbereiding op het vierkante geval.
Niet duidelijk hoe hij gedragen wordt.  Zeg hij ligt op rubber steunpunten waarvan de positie is aangepast.
Of je doet het in de ruimte, zonder zwaartekracht. Dan heeft het geen zin er zand op te strooien.
De algemene oplossing zonder randvoorwaarden is (gevonden op http://www.phy.davidson.edu/StuHome/timv/IntLab/ChlaRes/Theory/Theo.htm )

Clear[a,b,c,d,k]

u[x_,a_,b_,c_,d_,k_]=a Cosh[k x]+b Sinh[k x]+c Cos[k x]+d Sin[k x]

D[u[x,a,b,c,d,k],{x,4}]

Simplify[D[u[x,a,b,c,d,k],{x,4}]==k^4  u[x,a,b,c,d,k]]

Iets leesbaarder:

D[u[x],{x,4}]==k^4  u[x]

Dus vierde afgeleide van u is k^4  maal u. Dat is de differentiaalvergelijking. Zie bijvoorbeeld Rossing en Fletcher, Principles of Vibration and Sound.
Voor een plaat moeten we straks de vierde afgeleide vervangen door tweemaal de laplaciaan.

Laten we eerst de oneven oplossingen zoeken. (oneven funkties hebben grafiek met (0,0) als symmetriepunt.)

We gooien er nu drie van de vier randvoorwaarden in, te vinden bij zelfde bron.

threecond=Solve[(D[u[x,a,b,c,d,k],{x,3}]/.x->1)==0&&(D[u[x,a,b,c,d,k],{x,2}]/.x->-1)==0&&(D[u[x,a,b,c,d,k],{x,3}]/.x->-1)==0,{a,b,c}]//Simplify

Omdat Mathematica het antwoord heeft veranderd, zetten we het oude maar hier neer.

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])}}

Dan zoeken we uit wat nummer vier zegt :

con=D[u[x,a,b,c,d,k],{x,2}]/.x->1/.threecond

Simplify[con/(d (k^2))]

Deze transcendente funktie moet dus nul worden voor een dunne balk met twee vrije uiteinden. Dat geeft voorwaarden voor k.

trans[k_]=-2 Sin[k]+2 Cos[k] Tanh[k]

Plot[trans[k],{k,0,8}]

De eerste triloplossing heet u1. Eerst de bijbehorende waarde van k.

k=k/.FindRoot[trans[k]==0,{k,3}]

d=1;u1[x_]=u[x,a,b,c,d,k]/.threecond/.a->1//First

Clear[k,d];
Plot[u1[x],{x,-1,1},PlotRange->All]

De tweede oplossing heet u2.

k=k/.FindRoot[trans[k]==0,{k,6}]

d=1;u2[x_]=u[x,a,b,c,d,k]/.threecond/.a->1//First

Clear[k,d];
Plot[u2[x],{x,-1,1},PlotRange->All]

Plot[trans[k],{k,5,10}]

k=k/.FindRoot[trans[k]==0,{k,10}]

d=1;u3[x_]=u[x,a,b,c,d,k]/.threecond/.a->1//First

Clear[k,d];
Plot[u3[x],{x,-1,1},PlotRange->All]

Plot[trans[k],{k,8,13}]

k=k/.FindRoot[trans[k]==0,{k,13}]

d=1;u4[x_]=u[x,a,b,c,d,k]/.threecond/.a->1//First

Clear[k,d];
Plot[u4[x],{x,-1,1},PlotRange->All]

Plot[trans[k],{k,11,17}]

k=k/.FindRoot[trans[k]==0,{k,16}]

d=1;u5[x_]=u[x,a,b,c,d,k]/.threecond/.a->1//First

Clear[k,d];
Plot[u5[x],{x,-1,1},PlotRange->All]

Plot[Sign[trans[k]],{k,0,20}]

De eerste zestien waarden van k:

onevenk=Table[k/.FindRoot[trans[k]==0,{k,4.3+j Pi}],{j,0,15}]

3.926602312047919`

onevenk/Pi

Ze liggen op voorspelbare plek.

trans[k]

Plot[ Tanh[k],{k,0,4}]

Nu dan op jacht naar de even oplossingen.

threecondeven=Solve[(D[u[x,a,b,c,d,k],{x,3}]/.x->1)==0&&(D[u[x,a,b,c,d,k],{x,2}]/.x->-1)==0&&(D[u[x,a,b,c,d,k],{x,3}]/.x->-1)==0,{a,b,d}]//Simplify

(D[u[x,a,b,c,d,k],{x,2}]/.x->1/.threecondeven)/(c (k^2) )//Simplify

transeven[k_]=-2 (Cos[k]+Coth[k] Sin[k])

Plot[Sign[transeven[k]],{k,0,7}]

evenk=Table[k/.FindRoot[transeven[k]==0,{k,2.3+j Pi}],{j,0,15}]

evenk/Pi

Nu code die even en oneven geval combineert.

uopl[m_]:=(uopl[m]=
If[OddQ[m],Function[x,(u[x,a,b,c,d,k]/.threecond//First)/.d->1/.k->onevenk[[(m-1)/2]]//Chop//Evaluate],Function[x,(u[x,a,b,c,d,k]/.threecondeven//First)/.c->1/.k->evenk[[m/2]]//Chop//Evaluate]])/;m>1

uopl[0]=Function[x,1];uopl[1]=Function[x,x];

De frequentie van de trillende staaf is op een faktor na gegeven door

freqstaaf[0]=freqstaaf[1]=0

freqstaaf[m_]:=If[OddQ[m],(k/.k->onevenk[[(m-1)/2]])^2,(k/.k->evenk[[m/2]])^2]/;m>1

De faktor hangt onder meer van de stijfheid van de staaf af.

Vierkant

We gebruiken het voorgaande. Zie Mary Waller pag 122 voor de differentiaalverglijking.

Nu dan echt Chladni?

ContourPlot[u2[x]u3[y],{x,-1,1},{y,-1,1},Contours->{-0.01,0.01},PlotPoints->100]

Het is geen wonder dat een productoplossing er zo uit ziet. Zijn mijn randvoorwaarden wel goed? Niet echt waarschijnlijk. Een beetje naief: gewoon weer de tweede en derde afgeleide in de normale richting (naar buiten toe) nul stellen. Dat doen die produktoplossingen, dus daar houd ik het maar op. Volgens Lord Rayleigh, Theory of Sound, pagina 372 is het veel te moeilijk om het correct te doen. De plaatjes suggereren dat we er niet ver naast zitten. Technisch gezegd negeren we geloof ik de constante van Poisson, zie Rossing en Fletcher pag 76. Merk op dat u2[x]u3[y] niet eens aan de differentiaalvergelijking bij een vaste frequentie voldoet! Alle interactie tussen de x en de y-richting is verwaarloosd. Maar helemaal onzin is het idee niet: bij vaste x zie je de oplossing voor een staaf met vrije uiteinden, en bij vaste y ook weer.

ContourPlot[u1[x]u1[y],{x,-1,1},{y,-1,1},Contours->{-0.01,0.01},PlotPoints->100]

De frequentie van de plaat is , weer op een faktor na,

freqvierkant[m_,n_]=freqstaaf[m]+freqstaaf[n]

Althans, dat zou hij zijn als je alleen cosinussen en sinussen had. Dat wil zeggen, als de trilling werd gegeven door Cos[2 Pi freq t] Sin[k1 x] Sin[k2 y] 
bijvoorbeeld. Die trilling voldoet niet aan randvoorwaarden die bij een vrije rand horen, maar hij voldoet wel in het inwendige aan de vergelijking. En zijn
frequentie is de som van de bij Sin[k1 x] respectievelijk Sin[k2 y] behorende frequenties.

freqvierkant[2,2]

Union[Flatten[Table[freqvierkant[m,n]/freqvierkant[2,2],{m,6},{n,6}]]]

ListPlot[Union[Flatten[Table[freqvierkant[m,n]/freqvierkant[2,2],{m,6},{n,6}]]],PlotStyle->PointSize[0.02]]

Table[Plot[uopl[m][t],{t,-1,1},PlotLabel->"opl("<>ToString[m]<>") freq="<>ToString[freqstaaf[m]]],{m,0,10}]

ContourPlot[uopl[3][x]uopl[4][y],{x,-1,1},{y,-1,1},Contours->{-0.01,0.01},PlotPoints->100]

Deze mag niet van Mary Waller. Toch best symmetrisch.

ContourPlot[uopl[3][x]uopl[4][y]+uopl[4][x]uopl[3][y],{x,-1,1},{y,-1,1},Contours->{-0.01,0.01},PlotPoints->100,PlotLabel->"("<>ToString[3]<>","<>ToString[4]<>")plus("<>ToString[4]<>","<>ToString[3]<>") freq="<>ToString[freqvierkant[3,4]]]

Deze mag wel van Mary Waller.

ContourPlot[uopl[3][x]uopl[4][y]-uopl[4][x]uopl[3][y],{x,-1,1},{y,-1,1},Contours->{-0.01,0.01},PlotPoints->100]

Deze mag wel van Mary Waller, maar de volgende dus weer niet.

ContourPlot[uopl[3][x]uopl[4][y]-.5uopl[4][x]uopl[3][y],{x,-1,1},{y,-1,1},Contours->{-0.01,0.01},PlotPoints->100]

ContourPlot[uopl[7][x]uopl[4][y]-uopl[7][y]uopl[4][x],{x,-1,1},{y,-1,1},Contours->{-0.01,0.01},PlotPoints->100]

ContourPlot[uopl[9][x]uopl[4][y]-uopl[9][y]uopl[4][x],{x,-1,1},{y,-1,1},Contours->{-0.01,0.01},PlotPoints->100]

ContourPlot[uopl[11][x]uopl[4][y]-uopl[11][y]uopl[4][x],{x,-1,1},{y,-1,1},Contours->{-0.01,0.01},PlotPoints->100]

ContourPlot[Abs[uopl[11][x]uopl[4][y]-uopl[11][y]uopl[4][x]],{x,-1,1},{y,-1,1},Contours->{0.3},PlotPoints->100]

Mooie plaatjes toch? Deze randvoorwaarden zijn zo slecht nog niet.

Ronde schijf

We pogen nu iets aan het ronde geval te doen. Noem de uitwijking  f0[r,phi] (op tijdstip nul) . Een keer Lapaciaan toepassen geeft f1:

Clear[y,f0,la,ka]

f1[r_,phi_]=laplace[f0[r,phi]]=D[f0[r,phi],{r,2}]+D[f0[r,phi],{r,1}]/r+D[f0[r,phi],{phi,2}]/(r^2)

Nog eens Laplace:

f2[r_,phi_]=D[f1[r,phi],{r,2}]+D[f1[r,phi],{r,1}]/r+D[f1[r,phi],{phi,2}]/(r^2)//Simplify

Dat moet een constant veelvoud van f0[r,phi] wezen. En die constante, die ik la noem, is evenredig met het kwadraat van de frequentie. De operator respecteert fourierontwikkeling naar phi, dus we zoeken oplossingen van de vorm

f0[r_,phi_]=y[r]Cos[n phi]

f2[r_,phi_]=D[f1[r,phi],{r,2}]+D[f1[r,phi],{r,1}]/r+D[f1[r,phi],{phi,2}]/(r^2)//Simplify

f2[r,phi]/f0[r,phi]==la//Simplify

DSolve[(((n^2) (-4+(n^2)) y[r]+r ((1+2 (n^2)) y'[r]+r (-(1+2 (n^2)) y''[r]+r (2 Derivative[3][y][r]+r Derivative[4][y][r]))))/((r^4) y[r]))==la,y[r],r]

Mathematica doet maar wat. Die Gamma's zijn vaak oneindig. Die schalen we er uit.

y[r_,n_,a_,b_,c_,d_,la_]=(y[r]/.DSolve[(1/((r^4) y[r]))((n^2) (-4+(n^2)) y[r]+r ((1+2 (n^2)) y'[r]+r (-(1+2 (n^2)) y''[r]+r (2 Derivative[3][y][r]+r Derivative[4][y][r]))))==la,y[r],r])/.{C[1]-> a(((-1))^((2-n)/4))I^n/ Gamma[2-n],C[2]->b(((-1))^(-n/4))I^n/Gamma[1-n],C[3]->c(((-1))^(n/4))I^n/Gamma[1+n],C[4]->I^n d(((-1))^((2+n)/4))/Gamma[2+n]}//First

Dit kan eenvoudiger.

y[r_,n_,a_,b_,c_,d_,la_]=a BesselI[-n,(la^(1/4)) r]+b BesselJ[-n,(la^(1/4)) r]+c BesselI[n,(la^(1/4)) r]+d BesselJ[n,(la^(1/4)) r]

Als n=0 dan zijn dit er wat weinig? Hoe dan ook, we zoeken oplossingen die horen bij gladde oplossingen op de schijf. Dus ook glad als we tot een diagonaal beperken en bij het midden kijken. Dus laten we y[r]=r^mu + hogere orde invullen.Weglaten van hogere orde komt hier neer op weglaten van het rechterlid.

Clear[y];(((n^2) (-4+(n^2)) y[r]+r ((1+2 (n^2)) y'[r]+r (-(1+2 (n^2)) y''[r]+r (2 Derivative[3][y][r]+r Derivative[4][y][r]))))/( y[r]))/.{y[r]->r^mu,Derivative[m_][y][r]->D[r^mu,{r,m}]}//Expand

Nu de zogenaamde indiciaalvergelijking oplossen.

Solve[4 (mu^2)-4 (mu^3)+(mu^4)-4 (n^2)+4 mu (n^2)-2 (mu^2) (n^2)+(n^4)==0,mu]

Als n=0 dan is er dus alleen mu=0 en mu=2. Dus twee onafhankelijke oplossingen is genoeg. En als n=1 dan weer twee, omdat de negatieve afvalt. Als n=2 dan valt de oplossing met niet triviale constante term af, omdat je anders kreukels krijgt. en voor hogere n houd je er twee over die mu positief hebben. Eerst het geval n=0.

y[r_,0,a_,b_,la_]=a BesselI[0,(la^(1/4)) r]+b BesselJ[0,(la^(1/4)) r]

y[r_]=y[r,0,a,b,la]

y[r]+O[r]^6

f1[r,phi]/.{n->0,r->1}

We willen dat de Laplaciaan dubbel nul is bij straal 1. Daarmee bedoel ik dat ook zijn gradient er nul is. Dat lijkt nu een handig analogon van de eisen over de tweede en derde afgeleide voor het geval van de staaf. Zal ook wel niet echt de goede voorwaarde zijn. Maar de plaatjes worden mooi. En ik ben niet de eerste: Mathieu liep al tegen deze combinatie aan. 
Zie pagina 370 van Rayleigh. Hij is ook heel logisch als je in termen van "positief definiet symmetrisch" denkt.

y[r_]=y[r,0,a,b,la]/.Solve[(f1[r,phi]/.{n->0,r->1})==0,{a}]/.b->1//First

Nou willen we hem weer dubbel nul, dus stel ook de afgeleide nul op de rand Deze afgeleide wordt genomen in de richting "naar buiten". (Of stel de gradient nul.)

err[ka_,0]=D[f1[r,phi]/.n->0,r]/.r->1/.(la^(m_))->ka^(4m)

Plot[Sign[err[ka,0]],{ka,0,20}]

(ka/.Table[FindRoot[err[ka,0],{ka,n Pi}],{n,6}])/Pi

Dus heel dicht bij veelvouden van Pi. De bijbehorende y[r] zien er zo uit:

y0[r_]=y[r];

opl[0]=Table[la=(ka/.FindRoot[err[ka,0],{ka,n Pi}])^4;y0[r],{n,13}]

Table[Plot[opl[0][[i]]//Evaluate,{r,0,1}],{i,6}]

Bijbehorende Chladni patronen. Het grijze stuk hoort er niet bij. Dat is lucht.

Table[ContourPlot[(opl[0][[i]]/.r->Sqrt[u^2+v^2]//Evaluate ) If[u^2+v^2<1,1,0],{u,-1,1},{v,-1,1},Contours->{-0.01,0.01},PlotPoints->100],{i,6}]

Tot zover het geval n=0. Nu n>0.

Clear[la];y[r_,n_,a_,b_,c_,d_,la_]=a BesselI[-n,(la^(1/4)) r]+b BesselJ[-n,(la^(1/4)) r]+c BesselI[n,(la^(1/4)) r]+d BesselJ[n,(la^(1/4)) r]

y[r,1,a,b,c,d,la]+O[r]^6

Een oneven funktie waarvan we alleen het deel met r>0 verwerken.

Table[BesselI[n,r]==BesselI[-n,r]//Simplify,{n,6}]

Table[BesselJ[n,r]==BesselJ[-n,r]//Simplify,{n,6}]

Dus alweer zijn er minder dan vier gevonden. Mathematica doet maar wat. Voor ons is het genoeg.

y[r_]=y[r,nn,0,0,c,d,la]

f1[r,phi]/.{r->1}

We willen dat de Laplace dubbel nul is bij straal 1. Eerst maar eens nul.

y[r_]=y[r,nn,0,0,c,d,la]/.Solve[(f1[r,phi]/.{r->1,n->nn})==0,{c}]/.d->1//First

Nou dubbel nul. Dus afgeleide ook nog nul.

err[ka_,nn_]=D[f1[r,phi],r]/.r->1/.(la^(m_))->ka^(4m)/.phi->0/.n->nn

err[.7,1]

Plot[err[ka,1],{ka,0,4}]

Plot[Sign[err[ka,1]],{ka,0,20}]

(ka/.Table[FindRoot[err[ka,1],{ka,2+n Pi}],{n,6}])/Pi

Dus heel dicht bij Pi/2 plus veelvouden van Pi. De bijbehorende y[r] zien er zo uit.

opl[1]=Table[la=(ka/.FindRoot[err[ka,1],{ka,2+m Pi}])^4;(y[r]/.{nn->1,n->1}),{m,6}]

Table[Plot[opl[1][[i]]//Evaluate,{r,0,1}],{i,6}]

Bijbehorende Chladni patronen. Het grijze stuk hoort er niet bij. Dat is lucht.

Table[ContourPlot[(opl[1][[i]]/.r->Sqrt[u^2+v^2]//Evaluate )Cos[Arg[u+I v]] If[u^2+v^2<1,1,0],{u,-1,1},{v,-1,1},Contours->{-0.001,0.001},PlotPoints->100],{i,6}]

Nu net zo iets voor n=2.

Clear[la];(y[0]/.nn->2)

Plot[Sign[err[ka,2]],{ka,0,20}]

(ka/.Table[FindRoot[err[ka,2],{ka,2+n Pi}],{n,6}])

(ka/.Table[FindRoot[err[ka,2],{ka,2+n Pi}],{n,6}])/Pi

opl[2]=Table[la=(ka/.FindRoot[err[ka,2],{ka,2+m Pi}])^4;y[r]/.{nn->2,n->2},{m,6}]

Algemener werkt, voor niet al te grote waarden

freqschijf[m_,1]:=
(freqschijf[m,1]=
(ka/.FindRoot[err[ka,m],{ka,(.38 m +1)Pi}])^2)/;0<=m<11

freqschijf[m_,1]:=(freqschijf[m,1]=(((ka/.FindRoot[err[ka,m],{ka,1+Sqrt[freqschijf[m-1,1]]}]))^2))/;m>10

freqschijf[m_,n_]:=freqschijf[m,n]=(ka/.FindRoot[err[ka,m],{ka,Sqrt[freqschijf[m,n-1]]+Pi}])^2

ListPlot3D[Table[Sqrt[freqschijf[m,n]],{m,0,9},{n,9}],ViewPoint->{-2.396, -1.967, 1.356}]

ListPlot3D[Table[Sqrt[freqschijf[m,n]]-(.43m+1.02n)Pi,{m,0,9},{n,9}],ViewPoint->{-2.396, -1.967, 1.356}]

.43Pi

ListPlot3D[Table[Sqrt[freqschijf[m,n]],{m,0,40},{n,40}],ViewPoint->{-2.396, -1.967, 1.356},PlotRange->All]

ListPlot3D[Table[Sqrt[freqschijf[m,n]]-m-n Pi,{m,0,40},{n,40}],ViewPoint->{-2.396, -1.967, 1.356},PlotRange->All]

ListPlot3D[Table[Sqrt[freqschijf[m,n]]-(.4m+1.03n)Pi,{m,0,40},{n,40}],ViewPoint->{-2.396, -1.967, 1.356},PlotRange->All]

opl[d_Integer,m_Integer]:=(opl[d,m]=Function[x,(la=freqschijf[d,m]^2;
If[d==0,y0[x],y[x]/.nn->d])//Evaluate])/;d≥0&&m>0

Table[Plot[opl[2][[i]]//Evaluate,{r,0,1}],{i,6}]

Bijbehorende Chladni patronen. Het grijze stuk hoort er niet bij. Dat is lucht.

Table[ContourPlot[(opl[2][[i]]/.r->Sqrt[u^2+v^2]//Evaluate )Cos[2Arg[u+I v]] If[u^2+v^2<1,1,0],{u,-1,1},{v,-1,1},Contours->{-0.001,0.001},PlotPoints->100],{i,6}]

DensityPlot[(opl[2][[3]]/.r->Sqrt[u^2+v^2]//Evaluate )Cos[2Arg[u+I v]] If[u^2+v^2<1,1,0],{u,-1,1},{v,-1,1},Mesh->False,PlotPoints->100]

De frequenties zijn weer op een faktor na te bepalen.

ListPlot[Table[freqschijf[0,n],{n,5}],PlotStyle->PointSize[0.02]]

ListPlot[Table[freqschijf[1,n],{n,5}],PlotStyle->PointSize[0.02]]

ListPlot[Table[freqschijf[2,n],{n,5}],PlotStyle->PointSize[0.02]]

freqschijf[0,1]

ListPlot[{Sequence@@Table[{n,freqschijf[i,n]}/freqschijf[0,1],{i,0,1},{n,5}]},PlotStyle->PointSize[0.02]]

ListPlot[{Sequence@@Table[{n,freqschijf[i,n]/freqschijf[0,1]},{i,0,2},{n,5}]},PlotStyle->PointSize[0.02]]

ListPlot[{Sequence@@Table[{n,freqschijf[i,n](1.02)^2/(freqschijf[0,1](.43i+1.02n)^2)},{i,0,6},{n,8}]},AxesOrigin->{-1,1}]

En dan kan het instructief zijn om de Laplaciaan te plotten.

lap[i_,j_]:=(g0[r,phi]=opl[i,j][r]Cos[i phi];D[g0[r,phi],{r,2}]+D[g0[r,phi],{r,1}]/r+D[g0[r,phi],{phi,2}]/(r^2))

ContourPlot[(lap[2,3]/.{r->Sqrt[u^2+v^2],phi->Arg[u+I v]}//Evaluate ),{u,-1,1},{v,-1,1},PlotPoints->30]

Hij is mooi plat nul op de rand.

Table[ContourPlot[(lap[p,q]/.{r->Sqrt[u^2+v^2],phi->Arg[u+I v]}//Evaluate ),{u,-1,1},{v,-1,1},Contours->{-0.1,0.1},PlotPoints->60],{p,0,0},{q,3}]

Table[ContourPlot[(lap[p,q]/.{r->Sqrt[u^2+v^2],phi->Arg[u+I v]}//Evaluate ),{u,-1,1},{v,-1,1},PlotPoints->30],{p,1,3},{q,1,3}]

En dan kan het instructief zijn om de afgeleide naar buiten van de Laplaciaan te plotten.

dlap[i_,j_]:=(g0[r,phi]=opl[i,j][r]Cos[i phi];D[D[g0[r,phi],{r,2}]+D[g0[r,phi],{r,1}]/r+D[g0[r,phi],{phi,2}]/(r^2),r])

Table[ContourPlot[(dlap[p,q]/.{r->Sqrt[u^2+v^2],phi->Arg[u+I v]}//Evaluate ),{u,-1,1},{v,-1,1},PlotPoints->30],{p,1,3},{q,1,3}]

dlap[2,1]/.r->1

Table[dlap[p,q]/.{r->1},{p,1,3},{q,1,3}]

ContourPlot[Abs[opl[7,3][r]Cos[7 phi]/.{r->Sqrt[u^2+v^2],phi->Arg[u+I v]}]//Evaluate,{u,-1,1},{v,-1,1},PlotPoints->100,Contours->{0.1}]

En dan nu een controle met de Laplaciaan in Cartesische coordinaten.

f73=opl[7,3][Sqrt[x x+y y]]Cos[7 2ArcTan[(y/(x+(Sqrt[(x^2)+(y^2)])))]]

lhs=D[f73,x,x,x,x]+2D[f73,x,x,y,y]+D[f73,y,y,y,y];

Plot3D[(freqschijf[7,3]^2 f73-lhs),{x,.5,1},{y,-1,1},PlotRange->All,AxesEdge->{{-1,-1},{-1,1},{1,-1}}]

Tot slot nog wat betere code.

We willen dat de code robuuster is. Daarom maken we een Package.

(*:Name: Chladni` *)
(*:Author: Wilberd van der Kallen, 2006 *)
(*:Mathematica Version:5.2 *)
(*:Package Version:1.0 *)
(*:Sources:Rossing en Fletcher,Principles of Vibration and Sound. *)
(*Lord Rayleigh, Theory of Sound. *)
(*Mary Waller, Chladni Figures, A Study in Symmetry,1962 *)
BeginPackage["Chladni`"]

uopl::usage="uopl[m] gives a snapshot of an oscillating thin bending bar with free ends at -1 and 1."

freqstaaf::usage="freqstaaf[m] is proportional with the frequency associated with  uopl[m]."

freqschijf::usage="freqschijf[m,n] is proportional with the frequency associated with opl[m,n]."

freqvierkant::usage="freqvierkant[m,n] = freqstaaf[m] + freqstaaf[n]"

opl::usage="opl[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]."

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

uopl[0]=Function[x,1]

uopl[1]=Function[x,x]

uopl[m_Integer]:=(uopl[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])

freqvierkant[m_,n_]=freqstaaf[m]+freqstaaf[n]

freqstaaf[0]=0

freqstaaf[1]=0

freqstaaf[m_Integer]:=If[OddQ[m],onevk[(m-1)/2]^2,evk[m/2]^2]/;m>1


freqschijf[m_Integer,1]:=(freqschijf[m,1]=(ka/.FindRoot[err[ka,m],{ka,(38*m/100+1)*Pi}])^2)/;Inequality[0,LessEqual,m,Less,11]

freqschijf[m_Integer,1]:=(freqschijf[m,1]=(ka/.FindRoot[err[ka,m],{ka,1+Sqrt[freqschijf[m-1,1]]}])^2)/;m>10

freqschijf[m_Integer,2]:=(freqschijf[m,2]=(ka/.FindRoot[err[ka,m],{ka,1+Sqrt[freqschijf[m-1,2]]}])^2)/;m>5

freqschijf[m_Integer,n_Integer]:=(freqschijf[m,n]=(ka/.FindRoot[err[ka,m],{ka,Sqrt[freqschijf[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

opl[d_Integer,m_Integer]:=(opl[d,m]=Function[x,(la=freqschijf[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