finite difference method scheme

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

[aktham mansi]

The Drichlet boundary conditions are u(1; theta) = 1 and u(2; theta) = 0

[Irina Ayelen Jimenez]

The problem is the index j+0.5 in line 36. you should use the mean between j and j+1 of r. That is (r(j)+r(j+1))/2

[Torsten]

@aktham mansi

Maybe it's intersting to have the analytical solution for comparison:

r = linspace(1,2,41);
theta = linspace(0,2*pi,81);
theta = theta(1:end-1);
[r,theta] = meshgrid(r,theta);
uana = -log(r)/log(2) + 1;
surf(theta,r,uana)

[aktham mansi]

@Torsten i draw the analytical solution. can you draw the finite difference solution?

[Torsten]

Same commands as above. What's the problem ?

[aktham mansi]

@Torsten

i want to draw the finite difference solution and the analytical solution. ( can you add section to draw the finite difference solution?)

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

n

r = linspace(1,2,41);

theta = linspace(0,2*pi,81);

theta = theta(1:end);

[r,theta] = meshgrid(r,theta);

uana = -log(r)/log(2) + 1;

[x,y]=pol2cart(theta,r);

surface(x,y,uana);

colorbar;

[Torsten]

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;

Answer 1

[VBBV]

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

[VBBV]

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

contourf(u);
axis([1 40 1 13])
colorbar

Is this you are looking for?

[Kaid Noureddine]

Same answer as VBBV...

replace r(j+0.5) by (r(j)+r(j+1))/2

[aktham mansi]

@VBBV

it will look like that

not rectangular, but polar

[aktham mansi]

@Kaid Noureddine , thank you, could you please draw the domain for me?

the domain should look like this

.