# Wave Equation on an Annular Domain

This code solves a wave equation on an annular domain. Please note that this script defines functions at the end, which is only supported by MATLAB 2016b or later. The method is Chebyshev pseudospectral in r. The method is matrix based Fourier in theta. In polar coordinates the equation is written as:

## Contents

```
clear, close all;
```

## Initialization

Grid sizes.

Nx=96*2; Ny=64;

Chebyshev in yc (the computational coordinate in the radial direction).

[Dyc,yc]=cheb(Ny);

Fourier in xc (the computational coordinate for theta).

dxc = 2*pi/Nx; xc = -pi + (1:Nx)'*dxc; column = [0 .5*(-1).^(1:Nx-1).*cot((1:Nx-1)*dxc/2)]; Dxc = toeplitz(column,column([1 Nx:-1:2]));

Physical domain parameters and scaling for the physical grid.

Rin=1; % Inner radius. Rout=2; % Outer radius. dR=Rout-Rin; % Difference in radii. theta=xc; Dtheta=Dxc; r=Rin+(yc+1)*0.5*dR; [rr,thth]=meshgrid(r,theta);

For the operator on RHS.

over=1./r; Dr=Dyc*(2/dR); overr=1./rr; overrr2=(1./rr).^2; Dthetatheta=Dtheta*Dtheta; Drr=Dr*Dr;

Definitions for the Cartesian coordinates.

xx=rr.*cos(thth); yy=rr.*sin(thth);

Initial conditions.

un=sech(((xx-1.5).^2+yy.^2)/0.25).^16; up=un;

Parameters for plotting.

Ngx=1000; Ngy=200; rg=linspace(Rin,Rout,Ngy+1); thg=linspace(-pi,pi,Ngx+1); [rrg,ththg]=meshgrid(rg,thg); xxg=rrg.*cos(ththg); yyg=rrg.*sin(ththg);

Time stepping parameters.

% numsteps=20; numouts=500; dt=0.001; dt2=dt*dt; c=1; c2=c*c; % The above is the version for "live" movies. numsteps=250; numouts=9; dt=0.001; dt2=dt*dt; c=1; c2=c*c; t=0; for ii=1:numouts

ung=interp2(rr,thth,un,rrg,ththg,'spline'); % I use spline because it can handle the extrapolation as well. % I also make sure the grid wraps around so as to show the whole % annular domain.

%figure(ii) % For publishable version. figure(1) % For live movies. clf pcolor(xxg,yyg,ung),shading flat,caxis([-1 1]*0.5) % I caxis with a smaller range of values. % I do this so that later times are clearer. title(t) drawnow

```
for jj=1:numsteps
```

t=t+dt;

rhs=un*Drr'+overr.*(un*Dr')+overrr2.*(Dthetatheta*un); % Notice I can use fairly reasonably sized matrices. % This is because I keep the u variable as a matrix. % This means I need to get my left vs right multiplications % correct. uf=2*un-up+dt2*c2*rhs; uf(:,1)=0; uf(:,end)=0; up=un; un=uf;

```
end
```

```
end
```

And one last picture.

%figure(numouts+1) figure(2) clf ung=interp2(rr,thth,un,rrg,ththg,'spline'); pcolor(xxg,yyg,ung),shading flat,caxis([-1 1]*0.5) title(t),drawnow

## Chebyshev function

CHEB compute D = differentiation matrix, x = Chebyshev grid

function [D,x] = cheb(N) if N==0, D=0; x=1; return, end x = cos(pi*(0:N)/N)'; c = [2; ones(N-1,1); 2].*(-1).^(0:N)'; X = repmat(x,1,N+1); dX = X-X'; D = (c*(1./c)')./(dX+(eye(N+1))); % off-diagonal entries D = D - diag(sum(D')); % diagonal entries end