VERY URGENT::::: My code is giving me this error: Index in position 2 is invalid. Array indices must be positive integers or logical values. Can you please resolve this for me? It is very urgent and due in two days!!!

%####********Analytical Solutions *********#####
%#### This is simply plugging in the results from problem 1
%### Constants
sigmaa=0.005;
D=0.6;
L = (D/sigmaa)^0.5;
T=45;
So=10^8;

%### Slab Problem
x=linspace(0,45,10000);
U = (D/(2*L)-1/4)/(D/(2*L)+1/4);
C = (-So*T/sigmaa)*(T/4+D)/(-(D/(2*L)-1/4)^2/(D/(2*L)+1/4)*exp(-T/L)+(D/(2*L)+1/4)*exp(T/L));
A = U*C;
flux_slab=[];
for i = 1:x
    ans1 = A*exp(-i/L)+C*exp(i/L)+(So*i^2)/sigmaa;
    flux_slab = ans1;
end

%### Cylindrical Problem
x2=linspace(0,45/2,10000);
flux_cyl=[];
for i = 1:x2
    ans1 = -So/(2*sigmaa*(0.5*besseli(0,T/(2*L))+(D/L)*besseli(1,T/(2*L))))*besseli(0,i/L)+(So/sigmaa);
    flux_cyl = ans1;
end

%#### Sphere Problem
flux_sph=[];
for i = 1:x2
    flux_sph = (-So*sinh(i/L)/(i*4*sigmaa*(sinh(T/(2*L))/(2*T)+(2*D/L)*(T/2*cosh(T/(2*L))-L*sinh(T/(2*L)))/(T^2))))+So/sigmaa;
end

flux_comp_slab_10 = Diffusion('slab', 45, 10, 0.005, 0.6, 10^8);
flux_comp_slab_100 = Diffusion('slab', 45, 100, 0.005, 0.6, 10^8);
flux_comp_slab_1000 = Diffusion('slab', 45, 1000, 0.005, 0.6, 10^8);
flux_comp_slab_10000 = Diffusion('slab', 45, 10000, 0.005, 0.6, 10^8);
flux_comp_cyl_10 = Diffusion('cylinder', 45/2, 10, 0.005, 0.6, 10^8);
flux_comp_cyl_100 = Diffusion('cylinder', 45/2, 100, 0.005, 0.6, 10^8);
flux_comp_cyl_1000 = Diffusion('cylinder', 45/2, 1000, 0.005, 0.6, 10^8);
flux_comp_cyl_10000 = Diffusion('cylinder', 45/2, 10000, 0.005, 0.6, 10^8);
flux_comp_sph_10 = Diffusion('sphere', 45/2, 10, 0.005, 0.6, 10^8);
flux_comp_sph_100 = Diffusion('sphere', 45/2, 100, 0.005, 0.6, 10^8);
flux_comp_sph_1000 = Diffusion('sphere', 45/2, 1000, 0.005, 0.6, 10^8);
flux_comp_sph_10000 = Diffusion('sphere', 45/2, 10000, 0.005, 0.6, 10^8);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PLOTTING OPTIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x3=linspace(0,45/2,10);
x4=linspace(0,45/2,100);
x5=linspace(0,45/2,1000);
x6=linspace(0,45/2,10000);
x7=linspace(0,45/2,100000);

figure(1);
plot(x, flux_slab);
hold on
plot(x3*2, flux_comp_slab_10);
plot(x4*2, flux_comp_slab_100);
plot(x5*2, flux_comp_slab_1000);
plot(x6*2, flux_comp_slab_10000);
hold off

figure(2);
plot(x, flux_cyl);
hold on
plot(x3*2, flux_comp_cyl_10);
plot(x4*2, flux_comp_cyl_100);
plot(x5*2, flux_comp_cyl_1000);
plot(x6*2, flux_comp_cyl_10000);
hold off

figure(2);
plot(x, flux_sph);
hold on
plot(x3*2, flux_comp_sph_10);
plot(x4*2, flux_comp_sph_100);
plot(x5*2, flux_comp_sph_1000);
plot(x6*2, flux_comp_sph_10000);
hold off


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DEFINITIONS

function A = SOR(x, So)
S = zeros(length(x),1);     % BLANK MATRIX
for i = 1:(length(S)-1)
    S(i)=So*(x(i)^2);
end
A=S;
return
end


function B = SORSC(x, So)
S = zeros(length(x),1);     % BLANK MATRIX
for i = 1:(length(s)-1)
    S(i)=So;
end
B=S;
return
end


function C = Diffusion(geometry, Size, N, SigmaA, D, So)
x = linspace(0,Size,N)     ;       % Set the size to iterate and number of nodes
delta = (x(2)-x(1));                 % Distance between nodes
L = zeros(length(x), length(x));   % Start with a blank matrix


% GEOMETRY TYPE
if geometry=='slab'                 % Depending on geometry make tridiagonal terms
    S = SOR(x, So) / D;               % Recieve source term from source def
    for i = 1:(length(x)-1)                 % Iterate over tri diagonal
        L(i,i) = 2/(delta^2)+SigmaA/D;   % Diagonal Terms
        L(i,i+1) = -1/(delta^2);         % Above diagonal term
        L(i,i-1) = -1/(delta^2);         % Below diagonal term
    end
% Set boundary conditions at ends (Vacuum conditions)
L(0,0) = 0.5*(0.5-D/delta);
L(0,1) = D/(2*delta);
L(length(x)-1,length(x)-1) = 0.5*(0.5-D/delta);
L(length(x)-1, length(x)-2) = D/(2*delta);


elseif geometry=='cylinder'
    S = SORSC(x, So) / D ;  % Recieve source term from source def
    for i = 1:(length(x)-1)  % Iterate over tri diagonal
        L(i,i) = (x(i+1)/(x(i)*(delta^2))+1/(delta^2)+SigmaA/D); % Diagonal Terms
        L(i,i+1) = -x(i+1)/(x(i)*delta^2);                       % Above diagonal term
        L(i,i-1) = -1/(delta^2);                                 % Below diagonal term
    end
% Set boundary conditions at ends (Vacuum conditions 2/ symmetry)
L(0,0) = -1/delta;
L(0,1) = 1/delta;
L(length(x)-1,length(x)-1) = D/(2*delta);
L(length(x)-1, length(x)-2) = 0.5*(0.5-D/delta);


elseif geometry=='sphere'
    S = SORSC(x, So) / D;       % Recieve source term from source def
    for i = 1:(len(x)-1)  % Iterate over tri diagonal
        % Diagonal Terms
        L(i,i) = (x(i+1)^2)/((x(i)^2)*(delta^2)) + 1/(delta^2) + SigmaA/D;
        L(i,i+1) = -x(i+1)^2 / (x(i)^2 * delta^2); % Above diagonal term
        L(i,i-1) = -1/(delta^2); % Below diagonal term
    end
        % Set boundary conditions at ends (Vacuum conditions w/ symmetry)
L(0,0) = -1/delta;
L(0,1) = 1/delta;
L(length(x)-1,length(x)-1) = D/(2*delta);
L(length(x)-1, length(x)-2) = 0.5*(0.5-D/delta);

else
    sprintf("Not a valid geometry option");
end
flux = linsolve(L,S);
C=flux;
return
end

Answer 1

[Rik]

There are many issues with your question. The most important one is this: your question is not urgent.

After that, you do an element-wise comparison for the input char arrays, instead of using strcmp (or switch-case). The root of your problem is that you try to write the offset elements from the diagonal, but don't check if those are actually inside the matrix.

You also seem to think Matlab is 0-indexed. The first array element in Matlab is not retrieved with a 0, but with a 1. I tried to fix the explicit mentions, but you should check them.

Additionally, you had a typo in your SORSC function: a lower case s instead of upper case, and a typo where you used len instead of length.

The code below still doesn't work as intended, but I can't correct that, because your loops are unclear what they should do. You overwrite the result every iteration, but it doesn't even start, because 1:x doesn't return a valid array to iterate over.

%####********Analytical Solutions *********#####
%#### This is simply plugging in the results from problem 1
%### Constants
sigmaa=0.005;
D=0.6;
L = (D/sigmaa)^0.5;
T=45;
So=10^8;
%### Slab Problem
x=linspace(0,45,10000);
U = (D/(2*L)-1/4)/(D/(2*L)+1/4);
C = (-So*T/sigmaa)*(T/4+D)/(-(D/(2*L)-1/4)^2/(D/(2*L)+1/4)*exp(-T/L)+(D/(2*L)+1/4)*exp(T/L));
A = U*C;
flux_slab=[];
for i = 1:x
    ans1 = A*exp(-i/L)+C*exp(i/L)+(So*i^2)/sigmaa;
    flux_slab = ans1;
end
%### Cylindrical Problem
x2=linspace(0,45/2,10000);
flux_cyl=[];
for i = 1:x2
    ans1 = -So/(2*sigmaa*(0.5*besseli(0,T/(2*L))+(D/L)*besseli(1,T/(2*L))))*besseli(0,i/L)+(So/sigmaa);
    flux_cyl = ans1;
end
%#### Sphere Problem
flux_sph=[];
for i = 1:x2
    flux_sph = (-So*sinh(i/L)/(i*4*sigmaa*(sinh(T/(2*L))/(2*T)+(2*D/L)*(T/2*cosh(T/(2*L))-L*sinh(T/(2*L)))/(T^2))))+So/sigmaa;
end
flux_comp_slab_10 = Diffusion('slab', 45, 10, 0.005, 0.6, 10^8);
flux_comp_slab_100 = Diffusion('slab', 45, 100, 0.005, 0.6, 10^8);
flux_comp_slab_1000 = Diffusion('slab', 45, 1000, 0.005, 0.6, 10^8);
flux_comp_slab_10000 = Diffusion('slab', 45, 10000, 0.005, 0.6, 10^8);
flux_comp_cyl_10 = Diffusion('cylinder', 45/2, 10, 0.005, 0.6, 10^8);
flux_comp_cyl_100 = Diffusion('cylinder', 45/2, 100, 0.005, 0.6, 10^8);
flux_comp_cyl_1000 = Diffusion('cylinder', 45/2, 1000, 0.005, 0.6, 10^8);
flux_comp_cyl_10000 = Diffusion('cylinder', 45/2, 10000, 0.005, 0.6, 10^8);
flux_comp_sph_10 = Diffusion('sphere', 45/2, 10, 0.005, 0.6, 10^8);
flux_comp_sph_100 = Diffusion('sphere', 45/2, 100, 0.005, 0.6, 10^8);
flux_comp_sph_1000 = Diffusion('sphere', 45/2, 1000, 0.005, 0.6, 10^8);
flux_comp_sph_10000 = Diffusion('sphere', 45/2, 10000, 0.005, 0.6, 10^8);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PLOTTING OPTIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x3=linspace(0,45/2,10);
x4=linspace(0,45/2,100);
x5=linspace(0,45/2,1000);
x6=linspace(0,45/2,10000);
x7=linspace(0,45/2,100000);
figure(1);
plot(x, flux_slab);
hold on
plot(x3*2, flux_comp_slab_10);
plot(x4*2, flux_comp_slab_100);
plot(x5*2, flux_comp_slab_1000);
plot(x6*2, flux_comp_slab_10000);
hold off
figure(2);
plot(x, flux_cyl);
hold on
plot(x3*2, flux_comp_cyl_10);
plot(x4*2, flux_comp_cyl_100);
plot(x5*2, flux_comp_cyl_1000);
plot(x6*2, flux_comp_cyl_10000);
hold off
figure(2);
plot(x, flux_sph);
hold on
plot(x3*2, flux_comp_sph_10);
plot(x4*2, flux_comp_sph_100);
plot(x5*2, flux_comp_sph_1000);
plot(x6*2, flux_comp_sph_10000);
hold off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DEFINITIONS
function A = SOR(x, So)
S = zeros(length(x),1);     % BLANK MATRIX
for i = 1:(length(S)-1)
    S(i)=So*(x(i)^2);
end
A=S;
return
end
function B = SORSC(x, So)
S = zeros(length(x),1);     % BLANK MATRIX
for i = 1:(length(S)-1)
    S(i)=So;
end
B=S;
return
end
function C = Diffusion(geometry, Size, N, SigmaA, D, So)
x = linspace(0,Size,N)     ;       % Set the size to iterate and number of nodes
delta = (x(2)-x(1));                 % Distance between nodes
L = zeros(length(x), length(x));   % Start with a blank matrix
% GEOMETRY TYPE
switch geometry
    case 'slab'                 % Depending on geometry make tridiagonal terms
        S = SOR(x, So) / D;               % Recieve source term from source def
        for i = 1:(length(x)-1)                 % Iterate over tri diagonal
            L(i,i) = 2/(delta^2)+SigmaA/D;   % Diagonal Terms
            if (i+1)<size(L,2)
                L(i,i+1) = -1/(delta^2);         % Above diagonal term
            end
            if (i-1)>0
                L(i,i-1) = -1/(delta^2);         % Below diagonal term
            end
        end
        % Set boundary conditions at ends (Vacuum conditions)
        L(1,1) = 0.5*(0.5-D/delta);
        L(1,2) = D/(2*delta);
        L(length(x),length(x)) = 0.5*(0.5-D/delta);
        L(length(x), length(x)-1) = D/(2*delta);
    case 'cylinder'
        S = SORSC(x, So) / D ;  % Recieve source term from source def
        for i = 1:(length(x)-1)  % Iterate over tri diagonal
            L(i,i) = (x(i+1)/(x(i)*(delta^2))+1/(delta^2)+SigmaA/D); % Diagonal Terms
            if (i+1)<size(L,2)
                L(i,i+1) = -x(i+1)/(x(i)*delta^2);                       % Above diagonal term
            end
            if (i-1)>0
                L(i,i-1) = -1/(delta^2);                                 % Below diagonal term
            end
        end
        % Set boundary conditions at ends (Vacuum conditions 2/ symmetry)
        L(1,1) = -1/delta;
        L(1,2) = 1/delta;
        L(length(x),length(x)) = D/(2*delta);
        L(length(x), length(x)-1) = 0.5*(0.5-D/delta);
    case 'sphere'
        S = SORSC(x, So) / D;       % Recieve source term from source def
        for i = 1:(length(x)-1)  % Iterate over tri diagonal
            % Diagonal Terms
            L(i,i) = (x(i+1)^2)/((x(i)^2)*(delta^2)) + 1/(delta^2) + SigmaA/D;
            if (i+1)<size(L,2)
                L(i,i+1) = -x(i+1)^2 / (x(i)^2 * delta^2); % Above diagonal term
            end
            if (i-1)>0
                L(i,i-1) = -1/(delta^2); % Below diagonal term
            end
        end
        % Set boundary conditions at ends (Vacuum conditions w/ symmetry)
        L(1,1) = -1/delta;
        L(1,2) = 1/delta;
        L(length(x),length(x)) = D/(2*delta);
        L(length(x), length(x)-1) = 0.5*(0.5-D/delta);
    otherwise
        error('Not a valid geometry option');
end
flux = linsolve(L,S);
C=flux;
return
end

[Sanchit Sharma]

Thanks a lot for looking into this. I am very grateful. I am fairly new to programming hence I am having a few issues with the code. I have attached the code I am trying to convert in to matlab. I really appreciate this.

[Sanchit Sharma]

I want to insert an array of size x and update it for the final plotting!

[Rik]

I don't know any Python, but it looks like you mistranslated some things:

array=[]
for i in x:
    ans=[formula]
    array.append(ans)

%should be this in Matlab:
array=[];
for i=x
    array=[array ans];
end

The other point I mentioned is indeed the case: Python is 0-indexed, so all indices should be incremented by 1 when translating Python to Matlab.

Answer 2

[per isakson]

See

• Debug a MATLAB Program
• Examine Values While Debugging

[Sanchit Sharma]

I have tried everything. I am still having this issue!

[Rik]

You didn't. There are still m-lint warnings when you copy this code block to the Matlab editor. At the very least you should have no m-lint warnings left. Also, the debugger would have helped you see that i is 1, so i-0 is the zeroth position, which doesn't exist in Matlab, hence the error.