# How can I solve the wave equation for a circular rod in cylindrical coordinates using finite differences?

I have a problem with the stability of finite difference method for the wave equation in cylindrical coordinates.

the equation is: $$\frac{\partial^2 \omega_n}{\partial r^2}+\frac{1}{r}\frac{\partial \omega_n}{\partial r}+\frac{\partial^2\omega_n}{\partial z^2}-\frac{\omega_n}{r^2}=\frac{1}{c^2}\frac{\partial^2\omega_n}{\partial t^2}$$

and when I run it in Matlab, I get the following result:

Does anyone know where the problem is?

I put the code here so that if anyone knows where the problem is, help me out.

%domain
clear;
Lx=10;
Ly=10;
dx=0.1;
dy=dx;
nx=fix(Lx/dx);
ny=fix(Ly/dy);
x=linspace(0,Lx,nx);
y=linspace(0,Ly,ny);
%time
T=10;
%varables
wn=zeros(nx,ny);
wnm1=wn; %w at time n-1
wnp1=wn; %w at tome n+1
%parameters
CFL=0.5; %CFL=c.dt/dx
c=1;
dt=(CFL*dx/c);

%dt=1e-6;

%%initial conditions
%%time stepping loop
t=0;
while(t<T)
%reflecting boundary conditions
wn(:,[1 end])=0;
wn([1 end],:)=0;

%solution
t=t+dt;
wnm1=wn; wn=wnp1; %save current and previous arrays

%sourse
wn(50,50)=dt^2*20*sin(30*pi*t/20);
for i=2:nx-1 for j=2:ny-1

wnp1(i,j)=2*wn(i,j)-wnm1(i,j)...
+CFL^2*(((wn(i+1,j)-2*wn(i,j)+wn(i-1,j))/dx^2+(1/(2*i*dx^2))*(wn(i+1,j)-wn(i-1,j))+(wn(i,j+1)-2*wn(i,j)+wn(i,j-1))/dy^2-wn(i,j)/(i*dx)^2));

end,end
%chek convergence

%visualize at selected steps
clf;
subplot(2,1,1);
imagesc(x,y,wn'); colorbar;caxis([-0.02 0.02])
title(sprintf('t=%.2f',t));
subplot(2,1,2);
mesh(x,y,wn'); colorbar;caxis([-0.02 0.02])
axis([0 Lx 0 Ly -0.05 0.05]);
shg; pause(0.01);
end

• You can just export the image in MATLAB. Have you reduced CFL? Feb 11 '18 at 12:12
• CFL=0.5 and i reduced this Coefficient for several times,but the answer are still unstable @nicoguaro Feb 11 '18 at 12:50