The Drichlet boundary conditions are u(1; theta) = 1 and u(2; theta) = 0
finite difference method scheme
47 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
discretization with uniform (r*theta) (81 * 41), Implement Jacobi,
the discretization equation is
i tried this code:
error: Array indices must be positive integers or logical values.
someone help me please
% laplace equation - 2D - Jacobi Method - Cylindrical / Polar
% Coordinates
% Dirichlet BC conditions - Constant properties at boundaries
clc
clear all
%%%%%%%%%%%%%%% Inputs
r_in=1; % Inside Radius of polar coordinates, r_in, say 1 m
r_out = 2; % Outside Radins of polar coordinates, r_out, say 2 m
j_max = 40; % no. of sections divided between r_in and r_out eg 80, 160,320
dr = (r_out - r_in)/j_max; % section length, m
%nr = j_max+1; % total no. of radial points =81 or 161 or 321
% total angle = 2*pi
i_max= 80; % no. of angle steps eg 40, 80, 160
dtheta= 2*pi/i_max; % angle step, rad
%Ur_in=1; %BC1
%Ur_out=0; %BC2
r = 1:dr:2;
theta=0:dtheta:2*pi;
[r,theta]=meshgrid(r,theta);
%%%%%initialize solution array
u=zeros(j_max+1,i_max+1);%%%%81*41 matrix
u_0=zeros(j_max+1,i_max+1);
u(1,:)=u(1,:)+1;
u(2,:)=u(2,:)+0;
beta=dr^2/dtheta^2;
n=1;
k=0;
%%% j index for radius r and i index for phi%%%%
while k==0
u_0=u;
k=1;
for i=2:80
for j=2:40
r(j)=1+(j-1)*dr;
theta(i)=dtheta/2+(i-1)*dtheta;
u(i,j)=(r(j+0.5)*u_0(i,j+1)+r(j-0.5)*u_0(i,j-1)+beta*u_0(i+1,j)+beta*u_0(i-1,j))/(r(j+0.5)+r(j-0.5)+2*beta);
if abs(u(i,j)-u_o(i,j))>(10^-5)
k=0;
end
end
end
n=n+1;
end
7 commentaires
Torsten
le 8 Avr 2022
Modifié(e) : Torsten
le 8 Avr 2022
As I said: If nothing is wrong with your solution, the following should work:
r = linspace(1,2,41);
theta = linspace(0,2*pi,81);
[r,theta] = meshgrid(r,theta);
uana = -log(r)/log(2) + 1;
[x,y]=pol2cart(theta,r);
figure(1)
surface(x,y,uana);
figure(2)
surface(x,y,u)
colorbar;
Réponse acceptée
VBBV
le 8 Avr 2022
clc
clear all
%%%%%%%%%%%%%%% Inputs
r_in=1; % Inside Radius of polar coordinates, r_in, say 1 m
r_out = 2; % Outside Radins of polar coordinates, r_out, say 2 m
j_max = 40; % no. of sections divided between r_in and r_out eg 80, 160,320
dr = (r_out - r_in)/j_max; % section length, m
%nr = j_max+1; % total no. of radial points =81 or 161 or 321
% total angle = 2*pi
i_max= 80; % no. of angle steps eg 40, 80, 160
dtheta= 2*pi/i_max; % angle step, rad
%Ur_in=1; %BC1
%Ur_out=0; %BC2
r = 1:dr:2;
theta=0:dtheta:2*pi;
[r,theta]=meshgrid(r,theta);
%%%%%initialize solution array
u=zeros(i_max+1,j_max+1);%%%%81*41 matrix
u_0=zeros(i_max+1,j_max+1);
u(1,:)=u(1,:)+1;
u(2,:)=u(2,:)+0;
beta=dr^2/dtheta^2;
n=1;
k=0;
%%% j index for radius r and i index for phi%%%%
while k==0
u_0=u;
k=1;
for i=2:80
for j=2:40
r(j)=1+(j-1)*dr;
theta(i)=dtheta/2+(i-1)*dtheta;
u(i,j)=(((r(j)+r(j+1))/2)*u_0(i,j+1)+((r(j)+r(j-1))/2)*u_0(i,j-1)+beta*u_0(i+1,j)+beta*u_0(i-1,j))/(((r(j)+r(j+1))/2)+((r(j)+r(j-1))/2)+2*beta);
if abs(u(i,j)-u_0(i,j))>(10^-5)
k=0;
end
end
end
n=n+1;
end
u
4 commentaires
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Function Creation dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!