Solving partial differential finite difference

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!
Ithumbnail_806621601_252643.jpg

Réponses (1)

Irina Ayelen Jimenez
Irina Ayelen Jimenez le 8 Oct 2020

0 votes

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

Catégories

En savoir plus sur Partial Differential Equation Toolbox dans Centre d'aide et File Exchange

Produits

Version

R2013b

Community Treasure Hunt

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

Start Hunting!

Translated by