波动方程课件.ppt

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、反变换可得达朗贝尔公式：达朗贝尔公式： 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其他当存在阻尼时，弦振动的振幅会不断减小，定解问题为：当存在阻尼时，弦振动的振幅会不断减小，定解问题为： 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 两端自由的弦的振动两端自由的弦的振动(有阻尼、驱动力有阻尼、驱动力) ; 00, ; 00,; 0, ; 0, 0)0( ; sin/cos 22txutxutlutulxtlxAuuautxxtxxtt当存在阻尼时，弦振动的振幅会不断减小；而存在驱当存

15、在阻尼时，弦振动的振幅会不断减小；而存在驱动力时，振幅会逐渐增大；定解问题为：动力时，振幅会逐渐增大；定解问题为： ninininininininiiiiinInInnnininininininininiftutuuucuudtuxtuuuuuuuftuutuuuxtauuu21111112111111122111122211 2 211 0 , 0 ; , 2- 2%ex6091; (p172)两端两端自由自由弦振动的差分解弦振动的差分解(有有阻尼、驱动力阻尼、驱动力)clear;II=20; N=750; L=1; T0=30; a=1; lamda=0; %(阻尼系数阻尼系数)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; %(驱动力驱动力)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 两端固定弦振动问题之五（质量不均匀）两端固定弦振动问题之五（质量不均匀） 如果弦的密度可变，则弦中张力也非常数，因此波动方如果弦的密度可变，则弦中张力也非常数，因此波动方程的形式应为：程的形式应为：若密度和张力均为指数型：若密度和张力均为指数型：则将其带入后得到：则将其带入后得到： (假定初速为假定初速为0)方程的差分格式为：方程的差分格式为： 22xuxTxtux00 , ;/*0

19、, 220022xulxtriAxuxuxuTtut xxeTxTex 0 0 ;111111112020121112020112 / 212 / 2iiiiiiininininininininiuuxuuuxtTuuuuxuuuxtTuuu%ex611; (p176) 两端固定的弦波动的差分解两端固定的弦波动的差分解(质量密度不均匀质量密度不均匀)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、，边界固定，膜的初始形状为旋转抛物的圆膜，边界固定，膜的初始形状为旋转抛物面，而初始速度为面，而初始速度为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的匀质圆柱，高为的匀质圆柱，高为L，上、下底面固定，侧面自由，上、下底面固定，侧面自由，初始位移为零初始位移为零, 初始速度为初始速度为u02z, 求柱体内的振动情况。求柱体内的振动情况。定解问题为：

34、定解问题为：其解析解为：其解析解为：7.3 7.3 三维波动问题三维波动问题7.3.1 柱体内的振动柱体内的振动 ;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) 柱体内振动问题的差分解柱体内振动问题的差分解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 柱体外的振动问题之一（第一种零阶汉克尔函数）柱体外的振动问题之一（第一种零阶汉克尔函数） 圆柱面半径为圆柱面半径为0，其径向速度为，其径向速度为v=v0cos t, 求柱面在柱求柱面在柱外所辐射的声场

39、的速度势。外所辐射的声场的速度势。定解问题为：定解问题为：其解析解为：其解析解为： 0;0, 0;0,; cos, ;002tututvtuuauttt tieaHvitu 10002Re,%ex617;(p188) 柱面声源的解析解柱面声源的解析解(第一种第一种零阶汉克尔函数零阶汉克尔函数)；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 柱体外的振动问题之二（汉克尔函数的应用）柱体外的振动问题之二（汉克尔函数的应用） 圆柱面半径为圆柱面半径为0，其径向速度为，其径向速度为v=v0cos cos t, 求柱面求柱面在柱外所辐射的声场的速度势。在柱外所辐射的声场的速度势。定解问题为：定解问题为：其解析解为：其解析解为： 0;0, 0;0,; coscos, ;002tututvtuuauttt tieaHvaitu 11200 cos2Re,%ex618;(p190) 柱面声源的解析解柱面声源的解析解(第一种一阶汉克尔函数第一种一阶汉克尔函数)；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

43、,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的球面，径向速度分布为的球面，径向速度分布为v=v0coscost，试求，试求该球面所发射的恒稳声振动的速度势。该球面所发射的恒稳声振动的速度势。定解问题为：定解问题为：其解析解为：其解析解为：在远场近似下：在远场近似下： ; Recos, ; 002tirtevtrruuautarartarrrvearirrvtrutari

44、/ sin / coscos2 cos12Re,300 / 2300tararrvtru / sincos2,300%ex614;(p191) 偶极声源的解析解偶极声源的解析解；clear; k=20*3; r0=0.2; v0=2; a=2; w=k*a; theta=0:2*pi/50:2*pi; rho=0.2:0.1:4; Th,R=meshgrid(theta,rho); X,Y=pol2cart(Th,R); uu=0; dt=0.001;N=25;t=(0:N)*dt;for n=1:N+1; BB=cos(Th).*exp(i*(k*R-w*t(n); u=real(v0*r0

45、3/2*(-1./R.2+i*k./R).*BB); figure(1); subplot(1,2,1); surf(X,Y,u); view(-32,28); axis(-4,4,-4,4,-2.0,1.0); p(n)=getframe; subplot(1,2,2); contour(X,Y,u,21); axis(-4,4,-4,4); pause(0.25); q(n)=getframe;end;subplot(1,2,1); movie(p); subplot(1,2,2);movie(q);6.3.5 四极声源四极声源半径为半径为r0的球面，径向速度分布为的球面，径向速度分布为v

46、=v0(3cos 2+1)/4*cost，试求该球面所发射的恒稳声振动的速度势。试求该球面所发射的恒稳声振动的速度势。定解问题为：定解问题为：其解析解为：其解析解为： ; Recos, ; 2002tirtePvtrruuautararrtarPrrvePrarairrvitrutari / sin3/1 / coscos3 cos3/13Re,222400 / 2223400%ex616;(p196) 偶极声源的解析解；偶极声源的解析解；clear; k=60; r0=0.2; v0=2; a=2; theta=0:2*pi/50:2*pi; rho=0.2:0.1:4; Th,R=mesh

47、grid(theta,rho); X,Y=pol2cart(Th,R); P2=1/2.*(3*(cos(Th).2-1);for t=0:0.002:0.05; AA=i*v0*r04/3*(-1./R.3+i*k./R.2+k2/3./R); BB=P2.*exp(i*k*(r-a*t); u=real(AA.*BB); figure(1);surf(X,Y,u); view(-32,28); axis(-4,4,-4,4,-5,5); p(n)=getframe; figure(2);contour(X,Y,u,21); pause(0.2); axis(-4,4,-4,4);axis square;colorbar; q(n)=getframe; end; movie(p); movie(q);