Solving the 3D wave equation
This scripts solves the 3D wave equation using the FFT. It also illustrates how to create nice 3D graphics.
Contents
close all;clear;
Create the grid
N=64; z=linspace(-1,1,N+1); z=z(1:end-1); [xx,yy,zz]=meshgrid(z,z,z);
Parameters.
dt=0.001; L=2; dk=2*pi/L;
Wavenumbers
ksvec(1)=0; ksvec(N/2+1)=0; for ii=2:(N/2) ksvec(ii)=ii-1; ksvec(N/2+ii)=-N/2 + ii -1; end ksvec=ksvec*dk; [k,l,m]=meshgrid(ksvec,ksvec,ksvec); k2=k.*k; l2=l.*l; m2=m.*m;
Create the operators
lap=-(k2+l2+m2); myop=dt*dt*lap;
Create the ICs
r=sqrt(xx.*xx+yy.*yy+4*zz.*zz); u0 = exp(-((r-.4)/.1).^2); u=u0; um1=u0;
Time parameters
numsteps=100; numouts=10; t=0;
Plot ICs
figure(1) clf hold on p=patch(isosurface(xx,yy,zz,u,0.05)); set(p,'FaceColor','blue','EdgeColor','none','facealpha',.3) daspect([1 1 1]) view(3); camlight, lighting phong sl = slice(xx,yy,zz,u,0,[],[]); minc = max(min(u(:)),-1); maxc = min(max(u(:)),1); cbot = -abs(max(abs(minc),abs(maxc))); ctop = -cbot; caxis([cbot ctop]); colormap hot;grid on %colorbar set(sl,'edgecolor','none','facecolor','interp','ambientstrength',1,... 'diffusestrength',0,'specularstrength',0,'facealpha',0.6); mystr=['time = ' num2str(t)]; title(mystr); axis([-1 1 -1 1 -1 1]) drawnow

Main Loop
for ii=1:numouts
for jj=1:numsteps
Step time.
t=t+dt;
FFT in y.
uf=fft(u);
FFT in x.
uf=fft(uf,[],2);
FFT in z.
uf=fft(uf,[],3);
Make the laplacian piece.
uf=myop.*uf;
iFFT in y.
uf=ifft(uf);
iFFT in x.
uf=ifft(uf,[],2);
iFFT in z.
uf=ifft(uf,[],3);
Knock off the small imaginary parts.
uf=real(uf);
Take a leapfrog step forward.
up1=-um1+2*u+uf; um1=u; u=up1;
end
Step plot
figure(1) clf hold on p=patch(isosurface(xx,yy,zz,u,0.05)); set(p,'FaceColor','blue','EdgeColor','none','facealpha',.3) daspect([1 1 1]) view(3); camlight, lighting phong sl = slice(xx,yy,zz,u,0,[],[]); minc = max(min(u(:)),-1); maxc = min(max(u(:)),1); cbot = -abs(max(abs(minc),abs(maxc))); ctop = -cbot; caxis([cbot ctop]); colormap hot;grid on %colorbar set(sl,'edgecolor','none','facecolor','interp','ambientstrength',1,... 'diffusestrength',0,'specularstrength',0,'facealpha',0.6); mystr=['time = ' num2str(t)]; title(mystr); axis([-1 1 -1 1 -1 1]) drawnow










end