Solving partial differential finite difference
35 views (last 30 days)
Show older comments
Hi there,
I am trying to solve the wave equation in polar coordinates using finite difference method. The image below is the equation/finite difference. I tried to follow the method in here but I am not sure how to change x and y to theta and r on MATLAB. Any help would be greatly appriciated.
Thanks!
I
0 Comments
Answers (1)
Irina Ayelen Jimenez
on 8 Oct 2020
clear;clc;
RUN = 1;
% Domain
R = 2;
% number of elements
nx = 20; ny = 50;
dx = R/nx; dy = 2*pi/ny;
r = linspace(0.01,R,nx); phiy = linspace(0,2*pi,ny);
[Rx,phi] = meshgrid(r,phiy);
xx = r.*cos(phi); yy = r.*sin(phi);
%%% Variables
wn = zeros(nx,ny); % current
wnm1 = wn; % time n-1
wnp1 = wn; % time n+1
%%% Parameters
CFL = 0.5; % CFL = c.dt/dx
c = 1;
% dt = CFL*dx/c;
dt = dx*dy/(4*c);
% dt=0.01;
% choose steps
steps = 3000;
%%% Initial condition peak
% wn(3:6,2:7) = 1;
% initial condition 1st mode
m = 0;n = 1;
km = [0 1.842 3.052 4.2012;3.8317 5.33 6.70 8.0152; 7.0156 8.5363 9.9695 11.346];
kb=km(n+1,m+1);
p = cos(m*phi).*besselj(m,kb*Rx/R);
wn = p';
wnp1(:) = wn(:);
t = 0;
figure;(1)
tn = 1;
% wn = p; wnp1 = p;
if RUN == 1
for k = 1:steps
% Reflecting boundary condition
wn(end,:) = 0;
% wn(1,1) = 0;
% Solution
t = tn*dt;
wnm1 = wn; wn = wnp1; % save current and previus array
% Source
% wn(10,:) = dt^2*20*sin(10*pi*t/20);
% wn(:,ny) = wn(:,1);
for i = 2:nx-1, for j = 2:ny-1
ri = max([r(i),0.5*dx]);
wnp1(i,j) = c*dt^2*((wn(i+1,j)-2*wn(i,j)+wn(i-1,j))/dx^2 + 1/ri*(wn(i+1,j)-wn(i-1,j))/dx/2 + 1/ri^2*(wn(i,j+1)-2*wn(i,j)+wn(i,j-1))/dy^2) ...
+ 2*wn(i,j)-wnm1(i,j);
end, end
% for j = 2:ny-1
% i = nx;
% wnp1(i,j) = c*dt^2*(1/ri^2*(wn(i,j+1)-2*wn(i,j)+wn(i,j-1))/dy^2);
% end
wnp1(1,:) = wnp1(2,:); % center
wnp1(:,1) = wnp1(:,2); wnp1(:,ny) = wnp1(:,1); % angular continuity
clf;
surf(xx,yy,wn');title(sprintf('t = %.2f', t));
zlim([-2 2])
shg;pause(0.01)
Matrix(:,:,tn) = wn;
tn = tn + 1;
end
end
1 Comment
See Also
Categories
Find more on Visualization in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!