工程声学/汽车消声器：动画：代码
% 带扩大的管道 - 最后修改于 1999 年 11 月 3 日
% 由 Prof. L. Mongeau 允许转载
clear all;
clc;
close all;
Nsteps=30;
f=86*2; % 频率
rho=1.15;
c=344.;
rhoc=rho*c;
l=1; % 管道长度
a1=0.0762/2;
S1=pi*a1^2;
a2=0.1524/2;
S2=pi*a2^2;
m=S2/S1;
w=2*pi*f;
k=w/c;
t=linspace(2*pi/Nsteps/w,2*pi/w,Nsteps); % 每个周期 Nsteps 个时间步长
alpha1=2.93e-5/a1.*sqrt(f);
alpha2=2.93e-5/a2.*sqrt(f);
khat=k-i*alpha1;
kl=khat*l;
U0=0.01;
R=(1-m)/(m+1);
A=rho*c*U0./(exp(i*kl)-R*exp(-i*kl));
B=R*A;
C=A+B;
x=linspace(0,2,100); % 计算域
y=linspace(-0.2,0.2,32);
for m=1:50
for k=8:24
pt(k,m)=A*exp(i*khat*(l-x(m)))-B*exp(-i*khat*(l-x(m)));
end
end
for m=51:100
for k=1:32
pt(k,m)=C*exp(i*khat*(l-x(m)));
end
end
pmag=abs(pt);
pphase=angle(pt);
clim=max(max(pmag));
V=[-clim clim]; % 设置颜色方案
pt=pmag.*exp(i*pphase); % 红色：最大正值
% 蓝色：最小负压
figure(1)
pt2=real(pt);
H=pcolor(x,y,pt2);
axis equal % 消除失真
axis([0 2 -.2 .2])
M=moviein(Nsteps); % 用于 Matlab 电影动画
for k=1:Nsteps
pt2=real(pt*exp(i*w*t(k)));
H=pcolor(x,y,pt2);
axis equal
xlabel('x');
ylabel('y');
title('压力')
shading interp;
caxis(V);
axis([0 2 -.2 .2])
M(:,k)=getframe;
% 存储图像以供网络发布（位图）
% filnme=strcat('spher',int2str(k));
% eval(['print -dbitmap ' filnme ' -f']);
end
movie(M,20)