finite difference method scheme

75 views (last 30 days)
aktham mansi
aktham mansi on 8 Apr 2022
Edited: aktham mansi on 8 Apr 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 Comments
Torsten
Torsten on 8 Apr 2022
Edited: Torsten on 8 Apr 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;

Sign in to comment.

Accepted Answer

VBBV
VBBV on 8 Apr 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
u = 81×41
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 0 0.2246 0.3798 0.4896 0.5690 0.6277 0.6717 0.7054 0.7314 0.7518 0.7678 0.7804 0.7903 0.7981 0.8039 0.8082 0.8112 0.8129 0.8134 0.8129 0.8113 0.8087 0.8049 0.8000 0.7939 0.7863 0.7772 0.7663 0.7532 0.7377 0 0.1039 0.1974 0.2783 0.3468 0.4043 0.4521 0.4918 0.5246 0.5517 0.5739 0.5921 0.6068 0.6185 0.6276 0.6343 0.6389 0.6416 0.6424 0.6415 0.6390 0.6347 0.6287 0.6210 0.6115 0.5999 0.5863 0.5704 0.5520 0.5307 0 0.0637 0.1240 0.1796 0.2299 0.2747 0.3141 0.3486 0.3783 0.4038 0.4255 0.4437 0.4587 0.4709 0.4805 0.4877 0.4926 0.4955 0.4964 0.4953 0.4925 0.4877 0.4812 0.4729 0.4627 0.4506 0.4365 0.4205 0.4023 0.3820 0 0.0434 0.0850 0.1242 0.1606 0.1941 0.2244 0.2516 0.2757 0.2968 0.3152 0.3309 0.3441 0.3549 0.3635 0.3700 0.3745 0.3770 0.3778 0.3768 0.3740 0.3696 0.3636 0.3560 0.3467 0.3359 0.3235 0.3096 0.2941 0.2770 0 0.0309 0.0606 0.0889 0.1155 0.1403 0.1631 0.1838 0.2025 0.2191 0.2337 0.2463 0.2570 0.2658 0.2729 0.2783 0.2820 0.2841 0.2847 0.2837 0.2814 0.2777 0.2726 0.2662 0.2585 0.2496 0.2394 0.2282 0.2158 0.2023 0 0.0224 0.0441 0.0648 0.0843 0.1027 0.1196 0.1352 0.1494 0.1621 0.1733 0.1831 0.1914 0.1983 0.2039 0.2081 0.2110 0.2126 0.2130 0.2123 0.2104 0.2073 0.2033 0.1982 0.1921 0.1850 0.1771 0.1684 0.1588 0.1485 0 0.0165 0.0323 0.0476 0.0620 0.0756 0.0882 0.0998 0.1105 0.1200 0.1285 0.1359 0.1422 0.1475 0.1518 0.1550 0.1572 0.1584 0.1587 0.1581 0.1566 0.1542 0.1510 0.1471 0.1424 0.1370 0.1309 0.1243 0.1170 0.1092 0 0.0121 0.0238 0.0351 0.0457 0.0558 0.0651 0.0738 0.0817 0.0888 0.0952 0.1007 0.1055 0.1094 0.1126 0.1150 0.1167 0.1176 0.1178 0.1173 0.1161 0.1143 0.1119 0.1089 0.1054 0.1013 0.0967 0.0917 0.0863 0.0804 0 0.0089 0.0176 0.0259 0.0337 0.0412 0.0481 0.0545 0.0604 0.0657 0.0704 0.0745 0.0780 0.0810 0.0834 0.0851 0.0864 0.0870 0.0872 0.0868 0.0859 0.0845 0.0827 0.0804 0.0778 0.0747 0.0713 0.0676 0.0635 0.0592
  4 Comments
aktham mansi
aktham mansi on 8 Apr 2022
@Kaid Noureddine , thank you, could you please draw the domain for me?
the domain should look like this
.

Sign in to comment.

More Answers (0)

Categories

Find more on Graphics Objects in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by