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!
I
回答（1 个）
Irina Ayelen Jimenez
2020-10-8
0 个投票
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
类别
在 帮助中心 和 File Exchange 中查找有关 Partial Differential Equation Toolbox 的更多信息
2019-9-14
2022-4-8
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
