Why do I get the error "Subscript indices must either be real positive integers or logicals." ?

function GrCbvpthickn

% The effect of the solutal buoyancy parameter on the fluid velocity

% profiles is modelled.

clc;

%define space from 0 to 2

x=linspace(0,2,11);

GrC=0;

%initial guess for the answer

guess= bvpinit(x,[1;0;0;0;1;0;1;0]);

%setting options

options = bvpset('Stats','on');

hold on

for i = 1:5

GrC=GrC+1;

%calling bvp4c to solve the boundary value problem

sol=bvp4c(@bvp_func_g2,@boundary_conditions,guess,options,GrC);

%export solution

f = sol.y;

if GrC==1

f_1=f(1,:);

f_11=f(2,:);

X_1=sol.x;

else if GrC==2

f_2=f(1,:);

f_22=f(2,:);

X_2=sol.x;

else if GrC==3

f_3=f(1,:);

f_33=f(2,:);

X_3=sol.x;

else if GrC==4

f_4=f(1,:);

f_44=f(2,:);

X_4=sol.x;

else if GrC==5

f_5=f(1,:);

f_55=f(2,:);

X_5=sol.x;

end

end

end

end

end

end

l=length(f);

%plotting the solutions of f with GrT parameter

figure

plot(X_1,f_1,X_2,f_2,X_3,f_3,X_4,f_4,X_5,f_5)

legend('GrC=1','GrC=2','GrC=3','GrC=4','GrC=5')

title('The influence of solutal Grashof number on radial velocity profile')

xlabel('x')

ylabel('radial velocity,f(x)')

%plotting the solutions of f' with GrT parameter

figure

plot(X_1,f_11,X_2,f_22,X_3,f_33,X_4,f_44,X_5,f_55)

legend('GrC=1','GrC=2','GrC=3','GrC=4','GrC=5')

title('The influence of solutal Grashof number on axial velocity profile')

xlabel('x')

ylabel('axial velocity,df(x)/dx')

hold off

%----------------------------------------------------------------------------------------------------------

function dYdx = bvp_func_g2(x, Y, GrC)

%The system of first order ODEs is defined in this function with necessary

%values for parameters

%defining the constants

m=2;

J=1; %first second grade fluid parameter

L=1; %second second grade fluid parameter

M=4; %magnetic parameter

P=4; %porous parameter

GrT=2; %thermal grashof numbers

%GrC=1; %solatal grashof numbers

Pr=3; %prandtil number

Ec=1; %Eckert number

Sc=1; %Schimdit number

Kr=-0.3; %Chemical reaction paramter

alpha=1; %Heat generator or absorption paramter

H=1;

%defining derivatives

%Y(1)=f0, Y(2)=f1, Y(3)=f2,Y(4)=f3,Y(5)=theta, Y(6)=theta1,Y(7)=phi,

%Y(8)=phi1

dYdx(1)=Y(2);

dYdx(2)=Y(3);

dYdx(3)=Y(4);

dYdx(4)=1/((m-1)*x*Y(2)+((m+1)/2)*Y(1))*((1/J)*(m*Y(2)+((m-1)/sqrt(2))*x*Y(3)*Y(2)+((m+1)/2)*Y(3)*Y(1)+((m-1)/2)*x*Y(3)*Y(2)-(((m+1)^2)/2)*Y(4)*(L*Y(3))^(m/2)+M*Y(2)+P*Y(2)-GrT*Y(5)-GrC*Y(7))-(3*m-2)*Y(2)*Y(4)-((m-1)/2+((m-1)/sqrt(2))*x)*Y(3)*Y(4)-m*Y(4)*Y(1));

dYdx(5)=Y(6);

dYdx(6)=(2/(m+1))*(Pr*(Y(2)*Y(6)*((m-1)/2)*x+2*Y(5)*(((m+1)^2)/2)*(Y(1)+((m-1)/(m+1))*x*Y(1))-alpha*Y(5))-Pr*Ec(((m+1)/2)*Y(3)^2+J(m*Y(2)*Y(3)+x*((m-1)/2)*Y(3)^2+Y(2)*(Y(4)+Y(4)*x)+((m+1)/2)*(Y(2)*(1+((m-1)/(m+1)))+Y(3)*((m-1)/(m+1))))+H*(Y(1)*dYdx(4)+((m-1)/(m+1))*x*Y(2)*dYdx(4))+M*Y(2)^2));

dYdx(7)=Y(8);

dYdx(8)=(1/(m+1))*(Sc(Y(2)*Y(8)*x*(m-1)+(m+1)*(Y(1)+((m-1)/(m+1))*x*Y(2))+Kr*Y(7)));

dYdx=dYdx(:);

%-------------------------------------------------------------------------------------------------------------------------

function bc = boundary_conditions(Ya,Yb,GrC)

%boundary conditions at eta=1 and when eta goes to infinity is defined

%here as boundary a and b

%boundary conditions at x=0 and x=infinity

%boundary conditions at x=0

bc1=[Ya(1); Ya(2)-1; Ya(5)-1; Ya(7)-1];

%boundary conditions at x=infinity

bc2=[Yb(2); Yb(3); Yb(5); Yb(7)];

bc=[bc1;bc2];

this is my code and i get the error:

Subscript indices must either be real positive integers or logicals.

Error in GrCbvpthickn>bvp_func_g2 (line 107)

dYdx(6)=(2/(m+1))*(Pr*(Y(2)*Y(6)*((m-1)/2)*x+2*Y(5)*(((m+1)^2)/2)*(Y(1)+((m-1)/(m+1))*x*Y(1))-alpha*Y(5))-Pr*Ec(((m+1)/2)*Y(3)^2+J(m*Y(2)*Y(3)+x*((m-1)/2)*Y(3)^2+Y(2)*(Y(4)+Y(4)*x)+((m+1)/2)*(Y(2)*(1+((m-1)/(m+1)))+Y(3)*((m-1)/(m+1))))+H*(Y(1)*dYdx(4)+((m-1)/(m+1))*x*Y(2)*dYdx(4))+M*Y(2)^2));

Error in bvparguments (line 105)

testODE = ode(x1,y1,odeExtras{:});

Error in bvp4c (line 130)

bvparguments(solver_name,ode,bc,solinit,options,varargin);

Error in GrCbvpthickn (line 22)

sol=bvp4c(@bvp_func_g2,@boundary_conditions,guess,options,GrC);

Answer 1

[Image Analyst]

We get this asked every day. It's our most FAQQY of FAQs, so see the FAQ for a thorough discussion:

https://matlab.fandom.com/wiki/FAQ#%22Subscript_indices_must_either_be_real_positive_integers_or_logicals.%22

[Minoli Fernando]

I still don't know how to check whether all the indices are real positive integers in this problem. Can you please explain that

[Image Analyst]

EC starts out as a scalar, 1, but then you're trying to index it like it's a vector like

Ec(((m+1) ...............

was there supposed to be a * after Ec?

[Torsten]

Same for

J(m*Y(2)*Y(3) ...

[Minoli Fernando]

Thank you soo much. It cleared the error. But now i get the error singular jacobian encountered. What might be the course of that?

[Image Analyst]

Sorry, I have no idea about that error.

[Minoli Fernando]

Oh, that's okay your previous comment helped me. Thanks alot @Image Analyst @Torsten

Answer 2

[Star Strider]

Wioth the additional multiplicationo operators in ‘dYdx(6)’ and ‘dYdx(8)’ it appears to run without error —

%define space from 0 to 2
x=linspace(0,2,11);
GrC=0;
%initial guess for the answer
guess= bvpinit(x,[1;0;0;0;1;0;1;0]);
%setting options
options = bvpset('Stats','on');
% hold on 
for i = 1:5
    GrC=GrC+1;
   
    %calling bvp4c to solve the boundary value problem
    sol=bvp4c(@bvp_func_g2,@boundary_conditions,guess,options,GrC);
    %export solution
    f = sol.y;
    if GrC==1
        f_1=f(1,:);
        f_11=f(2,:);
        X_1=sol.x;
    else if GrC==2
        f_2=f(1,:);
        f_22=f(2,:); 
        X_2=sol.x;
    else if GrC==3
        f_3=f(1,:);
        f_33=f(2,:);
        X_3=sol.x;
    else if GrC==4
        f_4=f(1,:);
        f_44=f(2,:);
        X_4=sol.x;
    else if GrC==5
        f_5=f(1,:);
        f_55=f(2,:);
        X_5=sol.x;
        end
        end
        end
        end
    end
        
end

Warning: Unable to meet the tolerance without using more than 1250 mesh points.
The last mesh of 514 points and the solution are available in the output argument.
The maximum residual is 38.7022, while requested accuracy is 0.001.

The solution was obtained on a mesh of 1433 points.
The maximum residual is  3.870e+01. 
There were 58360 calls to the ODE function. 
There were 569 calls to the BC function. 

Warning: Unable to meet the tolerance without using more than 1250 mesh points.
The last mesh of 573 points and the solution are available in the output argument.
The maximum residual is 1261.26, while requested accuracy is 0.001.

Warning: Results may be inaccurate because of an ill-conditioned Jacobian whose reciprocal condition number is 2.74338e-17.

The solution was obtained on a mesh of 1702 points.
The maximum residual is  1.261e+03. 
There were 58507 calls to the ODE function. 
There were 571 calls to the BC function. 

Warning: Unable to meet the tolerance without using more than 1250 mesh points.
The last mesh of 800 points and the solution are available in the output argument.
The maximum residual is 57.341, while requested accuracy is 0.001.

Warning: Results may be inaccurate because of an ill-conditioned Jacobian whose reciprocal condition number is 2.86927e-32.

The solution was obtained on a mesh of 2062 points.
The maximum residual is  5.734e+01. 
There were 131026 calls to the ODE function. 
There were 910 calls to the BC function. 

Warning: Unable to meet the tolerance without using more than 1250 mesh points.
The last mesh of 795 points and the solution are available in the output argument.
The maximum residual is 624.938, while requested accuracy is 0.001.

Warning: Results may be inaccurate because of an ill-conditioned Jacobian whose reciprocal condition number is 1.82586e-20.

The solution was obtained on a mesh of 2211 points.
The maximum residual is  6.249e+02. 
There were 94304 calls to the ODE function. 
There were 657 calls to the BC function. 

Warning: Unable to meet the tolerance without using more than 1250 mesh points.
The last mesh of 1249 points and the solution are available in the output argument.
The maximum residual is 101.632, while requested accuracy is 0.001.

Warning: Results may be inaccurate because of an ill-conditioned Jacobian whose reciprocal condition number is 7.89545e-37.

The solution was obtained on a mesh of 3106 points.
The maximum residual is  1.016e+02. 
There were 237788 calls to the ODE function. 
There were 1022 calls to the BC function. 

    
    
    l=length(f);
%plotting the solutions of f with GrT parameter
figure
plot(X_1,f_1,X_2,f_2,X_3,f_3,X_4,f_4,X_5,f_5)
legend('GrC=1','GrC=2','GrC=3','GrC=4','GrC=5')
title('The influence of solutal Grashof number on radial velocity profile')
xlabel('x')
ylabel('radial velocity,f(x)')

%plotting the solutions of f' with GrT parameter
figure
plot(X_1,f_11,X_2,f_22,X_3,f_33,X_4,f_44,X_5,f_55)
legend('GrC=1','GrC=2','GrC=3','GrC=4','GrC=5')
title('The influence of solutal Grashof number on axial velocity profile')
xlabel('x')
ylabel('axial velocity,df(x)/dx')
hold off

%----------------------------------------------------------------------------------------------------------

function bc = boundary_conditions(Ya,Yb,GrC)
%boundary conditions at eta=1 and when eta goes to infinity is defined
%here as boundary a and b
%boundary conditions at x=0 and x=infinity
%boundary conditions at x=0
bc1=[Ya(1); Ya(2)-1; Ya(5)-1; Ya(7)-1];
%boundary conditions at x=infinity
bc2=[Yb(2); Yb(3); Yb(5); Yb(7)];
bc=[bc1;bc2];
end

function dYdx = bvp_func_g2(x, Y, GrC)
%The system of first order ODEs is defined in this function with necessary
%values for parameters
%defining the constants
m=2;
J=1; %first second grade fluid parameter
L=1; %second second grade fluid parameter
M=4; %magnetic parameter
P=4; %porous parameter
GrT=2; %thermal grashof numbers
%GrC=1; %solatal grashof numbers
Pr=3; %prandtil number
Ec=1; %Eckert number
Sc=1; %Schimdit number
Kr=-0.3; %Chemical reaction paramter
alpha=1; %Heat generator or absorption paramter
H=1;
%defining derivatives 
%Y(1)=f0, Y(2)=f1, Y(3)=f2,Y(4)=f3,Y(5)=theta, Y(6)=theta1,Y(7)=phi,
%Y(8)=phi1
dYdx(1)=Y(2);
dYdx(2)=Y(3);
dYdx(3)=Y(4);
dYdx(4)=1/((m-1)*x*Y(2)+((m+1)/2)*Y(1))*((1/J)*(m*Y(2)+((m-1)/sqrt(2))*x*Y(3)*Y(2)+((m+1)/2)*Y(3)*Y(1)+((m-1)/2)*x*Y(3)*Y(2)-(((m+1)^2)/2)*Y(4)*(L*Y(3))^(m/2)+M*Y(2)+P*Y(2)-GrT*Y(5)-GrC*Y(7))-(3*m-2)*Y(2)*Y(4)-((m-1)/2+((m-1)/sqrt(2))*x)*Y(3)*Y(4)-m*Y(4)*Y(1));
dYdx(5)=Y(6);
dYdx(6)=(2/(m+1))*(Pr*(Y(2)*Y(6)*((m-1)/2)*x+2*Y(5)*(((m+1)^2)/2)*(Y(1)+((m-1)/(m+1))*x*Y(1))-alpha*Y(5))-Pr*Ec*(((m+1)/2)*Y(3)^2+J*(m*Y(2)*Y(3)+x*((m-1)/2)*Y(3)^2+Y(2)*(Y(4)+Y(4)*x)+((m+1)/2)*(Y(2)*(1+((m-1)/(m+1)))+Y(3)*((m-1)/(m+1))))+H*(Y(1)*dYdx(4)+((m-1)/(m+1))*x*Y(2)*dYdx(4))+M*Y(2)^2));
dYdx(7)=Y(8);
dYdx(8)=(1/(m+1))*(Sc*(Y(2)*Y(8)*x*(m-1)+(m+1)*(Y(1)+((m-1)/(m+1))*x*Y(2))+Kr*Y(7)));
dYdx=dYdx(:);
end
%-------------------------------------------------------------------------------------------------------------------------

.