finite difference method scheme
62 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
aktham mansi
el 8 de Abr. de 2022
Editada: aktham mansi
el 8 de Abr. de 2022
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 comentarios
Torsten
el 8 de Abr. de 2022
Editada: Torsten
el 8 de Abr. de 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;
Respuesta aceptada
VBBV
el 8 de Abr. de 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 comentarios
Más respuestas (0)
Ver también
Categorías
Más información sobre Function Creation en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!