# drumhead problem, General case, Example 3 in Asmar 4.3 # u_{tt} = c^2(laplacian u), u=0 on clamped edge # circle radius=a, c=wave speed, # position = f(r,theta) = (a^2-r^2)*r*sin(theta), # velocity = g(r,theta) = (a^2-r^2)*r^2*sin(2*theta) a:=1;c:=1;N:=5; f:=(r,theta)->(a^2-r^2)*r*sin(theta); g:=(r,theta)->(a^2-r^2)*r^2*sin(2*theta); alpha1:=evalf([BesselJZeros(1,1..N)]); b1:=n->2*2*(a^5)*4/(alpha1[n]^3*BesselJ(2,alpha1[n])); u1:=(r,theta,t)->sin(theta)*sum(b1(n)*BesselJ(1,alpha1[n]*r)*cos(alpha1[n]*c*t),n=1..N); alpha2:=evalf([BesselJZeros(2,1..N)]); b2star:=n->24/(alpha2[n]^4*BesselJ(3,alpha2[n])); u2:=(r,theta,t)->sin(2*theta)*sum(b2star(n)*BesselJ(2,alpha2[n]*r)*sin(alpha2[n]*c*t),n=1..N); u:=(r,theta,t)->u1(r,theta,t)+u2(r,theta,t); addcoords(z_cylindrical,[z,r,theta],[r*cos(theta),r*sin(theta),z]); plot3d(f(r,theta),r=0..a,theta=0..2*Pi,coords=z_cylindrical,orientation=[13,81]); plot3d(g(r,theta),r=0..a,theta=0..2*Pi,coords=z_cylindrical); plots[animate](plot3d, [u(r,theta,t),r=0..a,theta=0..2*Pi,coords=z_cylindrical,orientation=[13,81]],t=0..3.1); # b1 can also be defined by integration directly. b2:=n->(2/(Pi*BesselJ(2,alpha1[n])^2))* int(int(f(r,theta)*BesselJ(1,alpha1[n]*r)*r*sin(theta),theta=0..2*Pi),r=0..a); seq(evalf(b2(j)-b1(j)),j=1..N);