# finite difference method scheme

53 views (last 30 days)
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.
% 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
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;

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
aktham mansi on 8 Apr 2022
@Kaid Noureddine , thank you, could you please draw the domain for me?
the domain should look like this
. R2018a

### Community Treasure Hunt

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

Start Hunting!