# 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
```