Solving partial differential finite difference

63 views (last 30 days)
Neyousha shahisavandi on 14 Sep 2019
Commented: aktham mansi on 8 Apr 2022
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

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
aktham mansi on 8 Apr 2022

Categories

Find more on PDE Solvers in Help Center and File Exchange

R2013b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by