1、7.1 一維波動問題一維波動問題7.2 二維波動問題二維波動問題7.3 三維波動問題三維波動問題trfuatu,222利用傅里葉變換可以得到：利用傅里葉變換可以得到：7.1 7.1 一維波動方程一維波動方程 7.1.1 無限長弦的自由振動無限長弦的自由振動 ;0, ;0,)( 2xtxuxtxuxuautxxtt atxatxdaatxatxtxu2121, ikdxexstepFikx10 ;0, ;0,sincos,22ktkUktkUaktkBaktkAtkUUikaUttt無限長弦的自由振動的定解問題為：無限長弦的自由振動的定解問題為：再用初始條件得：再用初始條件得：通過反變換可得通過

2、反變換可得達(dá)朗貝爾公式：達(dá)朗貝爾公式： akaktkaktktkU/sincos,一一. 解析解解析解1.若初始條件為：若初始條件為： 0 ; 07/47/3 ;/7sinxlxllxAx其他atxatxtxu21,則解析解則解析解為：為：%ex601; (p157) 一維無限長波動的解析解一維無限長波動的解析解(初速為初速為0)；clear; N=140; M=60; L=10; a=1; A=1; x1= 3*L/7;x2=4*L/7; x=L*(0:N)/N; dt=0.01; t=dt*(0:M); u0=A*sin(7*pi*x/L);u0(find(xx2)=0; figure(1

3、); h=plot(x,u0,linewidth,3);axis(0,L,-A,A); set(h,erasemode,xor); b=7*pi/L;for k=2:length(t); xl=x+a*t(k); xr=x-a*t(k); ul=A*sin(b*xl); ul(find(xlx2)=0; ur=A*sin(b*xr);ur(find(xrx2)=0; set(h,ydata,(ul+ur)/2 ); drawnow; pause(0.1);end;2. 初始條件為：初始條件為： 其他 ; 010 ; 10 xxx則解析解則解析解為：為： 1 ;2/1 10 ;2/0 ; 021

4、 1 ;2/1 10 ;2/0 ; 021212121,atxaatxaatxatxdaatxaatxaatxatxdadadadatxuatxatxatxatxatxatx其中：%ex602; (p159) 無限長弦波動的解析解無限長弦波動的解析解(初位移為初位移為0, 初速不為初速不為0)clear; M=100; N=80; a=1.0; L=10; T1=8; dx=L/M; dt=T1/N; x=-L:dx:L; t=0:dt:T1;X,T=meshgrid(x,t); xp=X+a*T; xp(find(xp=1)=1; xm=X-a*T; xm(find(xm=1)=1; S=(

5、xp-xm)/(2*a);figure(1); h=plot(x,S(1,:),linewidth,3); axis(-L L 0 .6); set(h,erasemode,xor); for k=2:N+1; pause(0.01); set(h,ydata,S(k,:); drawnow; end; ; 2-21211122211iiiiininininininixtuuxuuuuxtauuu二二. . 差分解差分解 ; 2 ;221211122211iiiiininininininixtuuxunuuuxtauuu%ex6012; (p157-161) 無窮長弦波動的差分解無窮長弦波動的

6、差分解clear; II=500; N=240; L=10; T=4; a=1.0; A=1; dx=L/II; dt=T/N; x=dx*(0:II); t=dt*(0:N); K=1:II+1;fai(K)=0; psi(K)=0; C=(a*dt/dx)2; I=2:II; %fai=A*sin(7*pi*x/L);fai(find(x4*L/7|x3*L/7)=0;% (初位移初位移) psi(find(x3*L/7)=1; %(初位移不為初位移不為0，初速為，初速為0)； u(1,:)=fai; u(2,:)=fai+dt*psi; figure(1); h=plot(x,u(1,:

7、),linewidth,3); %畫動畫畫動畫; axis(0,II*dx,-A,A); set(h,erasemode,xor); set(h,ydata,u(2,:); drawnow; pause(0.01); for k=2:N; u(k+1,1)=0; u(k+1,I+1)=0; u(k+1,I)=2*u(k,I)-u(k-1,I)+C*(u(k,I+1)-2*u(k,I)+u(k,I-1); set(h,ydata,u(k+1,:); drawnow; pause(0.01);end;figure(2); mesh(x,t,u);差分解程序差分解程序ex6012 差分解程序差分解程

8、序ex6012(無窮長，有限長無窮長，有限長) (無窮長，有限長無窮長，有限長)7.1.23 兩端固定的弦的自由振動兩端固定的弦的自由振動兩端固定的弦的自由振動的定解問題為：兩端固定的弦的自由振動的定解問題為： ;sin2;sin200lnlndxlxnxanBdxlxnxlA ;0, ;0,0, 0)0( 2xtxuxtxutlutulxuautxxtt解析解為：解析解為：其中其中1sinsincos,nnnlxnlatnBlatnAtxu1. 取取 l=1, a=1， ; 07/47/3 ;/7sin ; 0其他lxllxxx; 0 ; 7/17 ;/73sin7/74sin71 /73s

9、in7/74sin717nnBAn lnnnlnnnA2. 取取 l=1, a=1， ; 07/47/3 ; 1 ; 0其他lxlxx; 07/4cos7/3cos222nnAnnanlB%ex603; (p161) % 兩端固定的弦波動的解析解兩端固定的弦波動的解析解clear; N=100; M=500; K=100;a=1.0;L=1;T=.4;dx=L/M; dt=T/K; x=dx*(0:M); t=dt*(0:K); X,T=meshgrid(x,t); w=0;for n=1:N; p=(7+n)*pi/7; q=(7-n)*pi/7+eps; r=n*pi/L; s=n*pi;

10、A(n)=(sin(4*q)-sin(3*q)/q/7-(sin(4*p)-sin(3*p)/p/7; B(n)=0;%B(n)=2*a*L/(s*a)2*(cos(3*s/7)-cos(4*s/7);A(n)=0;w=w+(A(n)*cos(r*a*T)+B(n)*sin(r*a*T).*sin(r*X);end; ymax=1.1*max(max(abs(w);figure(1);h=plot(x,w(1,:),linewidth,3); axis(0 M*dx -ymax ymax); set(h,erasemode,xor);for k=2:K; pause(0.02); set(h,

11、ydata,w(k,:); drawnow; end;figure(2);mesh(X,T,w);解析解解析解 程序程序ex603 解析解解析解 程序程序ex603 7.1.4 兩端固定的弦的振動兩端固定的弦的振動(有阻尼有阻尼) ;00, )( ;07473 ;/7sin0,;0, ;0, 0)0( ;22xtxul/xl/lxxtxutlutulxuuauttxxtt其他當(dāng)存在阻尼時(shí)，弦振動的振幅會不斷減小，定解問題為：當(dāng)存在阻尼時(shí)，弦振動的振幅會不斷減小，定解問題為： nininiiinininininininiiiiiininininininininiuuucuututuuucuuux

12、tuuxuuutuuuxtauuu1112111110211111222112 *5 . 0 1/ 2 2 0; 2 ; 2- 2%ex6081; (p171)兩端固定弦振動的差分解兩端固定弦振動的差分解(有阻尼有阻尼)clear;II=500; N=300; L=10; T=4; a=1.0; A=1; lamda=5; dx=L/II; dt=T/N; x=dx*(0:II); t=dt*(0:N); K=1:II+1;fai(K)=0; psi(K)=0; C=(a*dt/dx)2; I=2:II; fai=A*sin(7*pi*x/L);fai(find(x4*L/7|x3*L/7)=

13、0; %(初位移初位移) u(1,:)=fai; u(2,I)=u(1,I)+0.5*C*(u(1,I+1)-2*u(1,I)+u(1,I-1);figure(1); h=plot(x,u(1,:),linewidth,3); %畫動畫畫動畫; axis(0,L,-1.1*A,1.1*A); set(h,erasemode,xor); set(h,ydata,u(2,:); drawnow; pause(0.01); for n=2:N; u(n+1,1)=0; u(n+1,II+1)=0; u(n+1,I)=2*u(n,I)-u(n-1,I)+C*(u(n,I+1)-2*u(n,I)+u(n

14、,I-1); u(n+1,I)=(u(n+1,I)+lamda*dt*u(n,I)/(1+lamda*dt); set(h,ydata,u(n+1,:); drawnow; pause(0.01);end;figure(2); mesh(x,t,u);title(lamda) =1/2 =1 =2 =5 =1/2 =1 =2 =5 7.1.5 兩端自由的弦的振動兩端自由的弦的振動(有阻尼、驅(qū)動力有阻尼、驅(qū)動力) ; 00, ; 00,; 0, ; 0, 0)0( ; sin/cos 22txutxutlutulxtlxAuuautxxtxxtt當(dāng)存在阻尼時(shí)，弦振動的振幅會不斷減?。欢嬖隍?qū)當(dāng)存

15、在阻尼時(shí)，弦振動的振幅會不斷減??；而存在驅(qū)動力時(shí)，振幅會逐漸增大；定解問題為：動力時(shí)，振幅會逐漸增大；定解問題為： ninininininininiiiiinInInnnininininininininiftutuuucuudtuxtuuuuuuuftuutuuuxtauuu21111112111111122111122211 2 211 0 , 0 ; , 2- 2%ex6091; (p172)兩端兩端自由自由弦振動的差分解弦振動的差分解(有有阻尼、驅(qū)動力阻尼、驅(qū)動力)clear;II=20; N=750; L=1; T0=30; a=1; lamda=0; %(阻尼系數(shù)阻尼系數(shù))dx=L/

16、II; dt=T0/N; x=dx*(0:II);t=dt*(0:N);X,T=meshgrid(x,t); K=1:II+1;fai(K)=0; psi(K)=0; f=zeros(N+1,II+1);u(1,:)=fai; u(2,:)=fai+dt*psi; C=(a*dt/dx)2; I=2:II; A=1; w=2; f=A*cos(pi*X/L).*sin(w*T)*dt2; %(驅(qū)動力驅(qū)動力)figure(1); h=plot(x,u(1,:),linewidth,3); %畫動畫畫動畫; axis(0,L,-.3,0.3); set(h,erasemode,xor); set(

17、h,ydata,u(2,:); drawnow; pause(0.01); for n=2:N; u(n+1,I)=2*u(n,I)-u(n-1,I)+C*(u(n,I+1)-2*u(n,I)+u(n,I-1); u(n+1,I)=(u(n+1,I)+lamda*dt*u(n,I)+f(n,I)/(1+lamda*dt); u(n+1,1)=u(n+1,2); u(n+1,II+1)=u(n+1,II); set(h,ydata,u(n+1,:); drawnow; pause(0.01);end;figure(2); mesh(x,t,u);兩端自由 ， f(x,t)=cos(pi*x)si

18、n(2*t) ， lamda=0兩端自由 ，f(x,t)=cos(pi*x)sin(2*t) ， lamda=2齊次邊界(兩端固定)，f(x,t)=cos(pi*x)sin(2*t) ，lamda=27.1.6 兩端固定弦振動問題之五（質(zhì)量不均勻）兩端固定弦振動問題之五（質(zhì)量不均勻） 如果弦的密度可變，則弦中張力也非常數(shù)，因此波動方如果弦的密度可變，則弦中張力也非常數(shù)，因此波動方程的形式應(yīng)為：程的形式應(yīng)為：若密度和張力均為指數(shù)型：若密度和張力均為指數(shù)型：則將其帶入后得到：則將其帶入后得到： (假定初速為假定初速為0)方程的差分格式為：方程的差分格式為： 22xuxTxtux00 , ;/*0

19、, 220022xulxtriAxuxuxuTtut xxeTxTex 0 0 ;111111112020121112020112 / 212 / 2iiiiiiininininininininiuuxuuuxtTuuuuxuuuxtTuuu%ex611; (p176) 兩端固定的弦波動的差分解兩端固定的弦波動的差分解(質(zhì)量密度不均勻質(zhì)量密度不均勻)clear; N=100; M=2000; L=1; T=1; a=10*sqrt(3); A=0.06; alfa=4;lamda=4; dx=L/N; dt=0.0005;x=dx*(0:N); t=dt*(0:M); C=(a*dt/dx)2

20、; I=2:N; J=1:N/2; d=2*A/L; u=zeros(N+1,M+1); K=N/2+1:N+1; u(J,1)=d*x(J); u(K,1)=d*(L-x(K); D=u(I+1,1)-2*u(I,1)+u(I-1,1); B=u(I+1,1)-u(I,1); u(I,2)=u(I,1)+C/2*(D+alfa*dx*B);figure(1); plot(0,1,0,0,r); hold on; %畫動畫畫動畫; h=plot(x,u(:,1),linewidth,3); axis(0,N*dx,-.1,.1); set(h,erasemode,xor);for k=2:M;

21、 set(h,ydata,u(:,k); drawnow; pause(0.01) D=u(I+1,k)-2*u(I,k)+u(I-1,k); B=u(I+1,k)-u(I,k); u(I,k+1)=2*u(I,k)-u(I,k-1)+C*(D+alfa*dx*B);end;7.1.7 7.1.7 非齊次邊界條件下弦的振動非齊次邊界條件下弦的振動; 0 ; 0 sin ; 0 ;2/20211111211kkknnKnnknknknknknkuuu;tAuuuuuxtauuu ; 00, ; 00,;sin, ; 0, 0)0( ;2txutxutAtlutulxuautxxtt 弦的一端弦的

22、一端(x=0)固定，另一端固定，另一端(x=l)作受迫振動作受迫振動Acos(t)，初，初位移和初速度均為位移和初速度均為0，其振動滿足下列定解問題：，其振動滿足下列定解問題：; 0 sin ; 0;2/ 5 . 0;2/ 21111111121211211knnKnkkkkknknknknknknku;tAuuuuuxtauuuuuxtauuu%ex6121; (p178) 非齊次邊界條件的弦波動的差分解非齊次邊界條件的弦波動的差分解clear; N=80; M=300; a=1; L=1; A=0.01; w=6; dx=L/N; dt=0.01; x=dx*(0:N); t=dt*(0:

23、M); C=(a*dt/dx)2; u=zeros(N+1,M+1); u(N+1,1:M+1)=A*sin(w*t); %初位移初位移； I=2:N; u(I,2)=u(I,1)+C/2*(u(I+1,1)-2*u(I,1)+u(I-1,1); figure(1); plot(0,1,0,0,r); hold on; %畫動畫畫動畫; h=plot(x,u(:,1),linewidth,3); axis(0,N*dx,-.1,.1); set(h,erasemode,xor);for k=2:M; set(h,ydata,u(:,k); drawnow; pause(0.01); u(I,k

24、+1)=2*u(I,k)-u(I,k-1)+C*(u(I+1,k)-2*u(I,k)+u(I-1,k);end;figure(2);mesh(x,t,u);差分解 解析解7.1.8 桿的縱振動桿的縱振動長為長為L的桿，可以沿縱向的桿，可以沿縱向x振動，初始位移振動，初始位移c0 x/L，初速為，初速為0，研究其振動。其定解問題是：研究其振動。其定解問題是：; 0 0 , ;/0 ,0, ; 0, 00,0 ;02xuLxcxu; tlututlxuautxxxxtt10/cos/cos,nnlxnlatnAAtxu1. 解析解解析解其中：其中：取取 c0 =0.05, a=1, l=1, N=

25、30,.3 , 2 , 1 ;124 ;222000nncAcAn%ex613; (p181) 桿的縱振動的解析解桿的縱振動的解析解；clear; N=50; M=100; K=1000; L=1;T0=5; a=1.0; C0=0.1; dx=L/M; dt=T0/M; x=dx*(0:M);t=dt*(0:M); X,T=meshgrid(x,t); w=1/2;for n=1:2:N; p=n*pi/L;w=w-4*L2/p2*cos(p*a*T).*cos(p*X);end; w=w*C0;figure(1); subplot(2,1,1); h1=plot(x,w(1,:),line

26、width,3); set(h1,erasemode,xor); axis(0,M*dx,0,C0); %畫動畫畫動畫; subplot(2,1,2); xx=1:5:M+1; yy=0*xx; h2=plot(x(xx),yy,r.,marker,.,markersize,25); set(h2,erasemode,xor); axis(0,1+C0,-C0,C0);for k=2:K+1; set(h1,ydata,w(k,:);drawnow; pause(0.2); set(h2,xdata,x(xx)+w(k,xx);drawnow; pause(0.2);end;7.2 7.2 二

27、維波動問題二維波動問題; 00 , ;sin0 ,0, ; 0, 0 , 0, ; 0, 00,0 ,0 ;2yxuy/baxAxyxu; tbxutxu;tyautyutbyaxuucutyyxxtt ;/sin ; 0 ;4/*5 . 0 ;4/21,1, 11 , 111,1, 11,11,1, 121,2,1, 1,1, 121,1,byaxAxuuuuuuuuuuxtcuuuuuuuxtcuuujiijinJinjIninjjijijijijijijinjinjinjinjinjinjinjinji7.2.1 矩形膜的振動矩形膜的振動 矩形膜的四周邊界固定，膜的初位移已知。定解問題是

28、：矩形膜的四周邊界固定，膜的初位移已知。定解問題是：顯式差分格式顯式差分格式：(取取 )1 2 1 ; 1 A;c;ba%ex621; (p183) 二維矩域波動方程的差分解二維矩域波動方程的差分解clear; N=80; M=40; K=500; a=2; b=1; c=1; A=1; dx=a/N; dy=b/M; dt=0.01; t=dt*(0:K); C=(c*dt/dx)2; x=dx*(0:N); y=dx*(0:M); Y,X=meshgrid(y,x); u1=zeros(N+1,M+1); u2=u1; u3=u1; I=2:N; J=2:M; u1=A*X.*(X-a).

29、*sin(pi*Y/b); %初位移初位移； u2(I,J)=u1(I+1,J)+u1(I,J+1)-4*u1(I,J)+u1(I-1,J)+u1(I,J-1); u2(I,J)=u1(I,J)+C/2*u2(I,J);figure(1);mesh(X,Y,u1); title(二維波動二維波動);axis(0 2 0 1 -1 1);p(1)=getframe; for k=2:K+1; u3(I,J)=u2(I+1,J)+u2(I,J+1)-4*u2(I,J)+u2(I-1,J)+u2(I,J-1); u3(I,J)=2*u2(I,J)-u1(I,J)+C*u3(I,J); u1=u2;

30、u2=u3;if mod(k,5)=0; mesh(X,Y,u3);axis(0 2 0 1 -1 1); p(k/5+1)=getframe; end;end; movie(p); ; 00 , ;10 , ; 0,0,0 ;/ 2000022tttu/uu tutuucucu; ;/1 ; 0; ; ; ) 1(2/2/ *5 . 0 ; ) 1(2/2/ 21220011112111111111112121111211iiiinInniiiiiiininininininininiuuuuuuuiuuuuutcuuiuuuuutcuuu7.2.2 圓膜的振動圓膜的振動 半徑為半徑為0的圓膜

31、，邊界固定，膜的初始形狀為旋轉(zhuǎn)拋物的圓膜，邊界固定，膜的初始形狀為旋轉(zhuǎn)拋物面，而初始速度為面，而初始速度為0，求膜的振動。定解問題是：，求膜的振動。定解問題是：顯式差分格式顯式差分格式：(取取 )4 . 0 1 ; 1 00u;a%ex622; (p186) 二維圓域波動方程的差分解二維圓域波動方程的差分解clear; N=90; K=500; a=1; r0=1; u0=0.4; dr=r0/N; dt=0.01; t=dt*(0:K); C=(a*dt/dr)2; r=dr*(0:N); fai=2*pi*(0:N)/N; R,Fai=meshgrid(r,fai);X,Y=pol2car

32、t(Fai,R); u1=zeros(1,N+1); u2=u1; u3=u1; I=2:N; u1=u0*(1-(r/r0).2); %初位移初位移； u2(I)=u1(I+1)-2*u1(I)+u1(I-1)+(u1(I+1)-u1(I-1)/(2*(I-1); u2(I)=u1(I)+C/2*u2(I); u2(1)=u2(2);figure(1);mesh(X,Y,ones(N+1,1)*u1); title(二維波動二維波動);axis(-r0 r0 -r0 r0 -.5 .5);p(1)=getframe; for k=2:K+1; u3(I)=u2(I+1)-2*u2(I)+u2

33、(I-1)+(u2(I+1)-u2(I-1)/(2*(I-1); u3(I)=2*u2(I)-u1(I)+C*u3(I); u3(1)=u3(2); u1=u2; u2=u3;if mod(k,5)=0; mesh(X,Y,ones(N+1,1)*u3); axis(-r0 r0 -r0 r0 -.5 .5);p(k/5+1)=getframe;end;end; movie(p); 半徑為半徑為0的勻質(zhì)圓柱，高為的勻質(zhì)圓柱，高為L，上、下底面固定，側(cè)面自由，上、下底面固定，側(cè)面自由，初始位移為零初始位移為零, 初始速度為初始速度為u02z, 求柱體內(nèi)的振動情況。求柱體內(nèi)的振動情況。定解問題為：

34、定解問題為：其解析解為：其解析解為：7.3 7.3 三維波動問題三維波動問題7.3.1 柱體內(nèi)的振動柱體內(nèi)的振動 ;0, ; 00,; 0, ; 0, 0, ; 0, ;2002zutzutzutLzutzutzuuauttt 22011101010212001/ 1/sinsin1 81,LmxkxJz/LmkatkaxJxmLutzunnmnnnm其中： ;0 , ; 00 , ; 0, ; 0, 0 , ; 0,0;0 ,0 ;/ 2000022zuzuzu tLututzutLzuuuauautzztt ; 0 ; 0; ; 0 ; 0; ; ; ) 1( 2/4* / *5 . 0

35、; ) 1( 2/4* / 21,11,11 ,1, 11, 21, 1111111,11111,1, 12201,2,111, 1,1, 121,1,jinJininjInjnjiijiiijijijijijinininjinjinjinjinjinjinjinjiuuuuuuiuuuuuuutatzuuuiuuuuuuutauuu顯式差分格式顯式差分格式：(取取 )4 . 0 1 ; 1 00u;a%ex615; (p187) 柱體內(nèi)振動問題的差分解柱體內(nèi)振動問題的差分解clear; N=50; M=150; K=1000; a=0.5; b=1.5; c=1; A=1; dr=a/N;

36、dz=b/M; dt=0.005; t=dt*(0:K); C=(c*dt/dz)2; r=dr*(0:N); z=dz*(0:M); u1=zeros(N+1,M+1); u2=u1; u3=u1; u0=A*(r).2*z; I=2:N; J=2:M; I1=1:N+1; J1=1:M+1; u2(I,J)=u1(I+1,J)+u1(I,J+1)-4*u1(I,J)+u1(I-1,J)+u1(I,J-1); dudr=(u1(I+1,J)-u1(I-1,J)./(2*(I-1)*ones(1,M-1); u2(I,J)=u1(I,J)+u0(I,J)*dt+C/2*(u2(I,J)+dud

37、r); u2(1,J1)=u2(2,J1);figure(1);mesh(z,r,u1);title(二維波動二維波動);axis(0 1.5 0 .5 -.04 .04); p(1)=getframe; mesh(z,r,u2);title(二維波動二維波動);p(2)=getframe;for k=2:K; u3(I,J)=u2(I+1,J)+u2(I,J+1)-4*u2(I,J)+u2(I-1,J)+u2(I,J-1); dudr=(u2(I+1,J)-u2(I-1,J)./(2*(I-1)*ones(1,M-1); u3(I,J)=2*u2(I,J)-u1(I,J)+C*(u3(I,J

38、)+dudr); u3(1,J1)=u3(2,J1); u3(N+1,J1)=0; u3(I1,1)=0; u3(I1,M+1)=0; u1(I,J)=u2(I,J); u2(I,J)=u3(I,J);if mod(k,5)=0; mesh(z,r,u3); axis(0 1.5 0 .5 -.04 .04); p(k/5+1)=getframe; end;end; movie(p);6.3.2 柱體外的振動問題之一（第一種零階漢克爾函數(shù)）柱體外的振動問題之一（第一種零階漢克爾函數(shù)） 圓柱面半徑為圓柱面半徑為0，其徑向速度為，其徑向速度為v=v0cos t, 求柱面在柱求柱面在柱外所輻射的聲場

39、的速度勢。外所輻射的聲場的速度勢。定解問題為：定解問題為：其解析解為：其解析解為： 0;0, 0;0,; cos, ;002tututvtuuauttt tieaHvitu 10002Re,%ex617;(p188) 柱面聲源的解析解柱面聲源的解析解(第一種第一種零階漢克爾函數(shù)零階漢克爾函數(shù))；clear; k=60; r1=1; r2=15; v0=1; a=2; w=6; theta=2*pi*(0:1/50:1); rho=r1:0.2:r2; Th,R=meshgrid(theta,rho); X,Y=pol2cart(Th,R); H=besselh(0,1,R); dt=0.05;

40、 N=40; t=(0:N)*dt;for n=1:N+1; u=real(-i*pi*v0*r1/2*H*exp(-i*w*t(n); figure(1);subplot(1,2,1);surf(X,Y,u); view(-32,28); axis(-12,12,-12,12,-2,2); p(n)=getframe; subplot(1,2,2);contour(X,Y,u,11); pause(0.2); axis(-15,15,-15,15); q(n)=getframe; end; subplot(1,2,1);movie(p); subplot(1,2,2);movie(q);6.

41、3.3 柱體外的振動問題之二（漢克爾函數(shù)的應(yīng)用）柱體外的振動問題之二（漢克爾函數(shù)的應(yīng)用） 圓柱面半徑為圓柱面半徑為0，其徑向速度為，其徑向速度為v=v0cos cos t, 求柱面求柱面在柱外所輻射的聲場的速度勢。在柱外所輻射的聲場的速度勢。定解問題為：定解問題為：其解析解為：其解析解為： 0;0, 0;0,; coscos, ;002tututvtuuauttt tieaHvaitu 11200 cos2Re,%ex618;(p190) 柱面聲源的解析解柱面聲源的解析解(第一種一階漢克爾函數(shù)第一種一階漢克爾函數(shù))；clear; k=60; r1=1; r2=15; v0=1; a=2; w=

42、6; theta=2*pi*(0:1/50:1); rho=r1:0.2:r2; Th,R=meshgrid(theta,rho); X,Y=pol2cart(Th,R); H=besselh(1,1,R).*cos(Th); dt=0.1; N=25; t=(0:N)*dt;for n=1:N+1; u=real(-i*pi*v0*r12*w/(2*a)*H*exp(-i*w*t(n); figure(1);subplot(1,2,1);surf(X,Y,u); axis(-12,12,-12,12,-2,2); p(n)=getframe; subplot(1,2,2);contour(X,Y,u,15); pause(0.1); axis(-15,15,-15,15); q(n)=getframe; end; subplot(1,2,1);movie(p); subplot(1,2,2);movie(q);6.3.4 偶極聲源偶極聲源 半徑為半徑為r0的球面，徑向