What is the meaning of "Index in position 1 exceeds array bounds (must not exceed 5)".

Hi everyone, I'm trying to run a code below:

syms a(x) b(x) c(x) d(x) e(x) A B C D E FF G H 
eqn1 = A * diff(a(x),x,2) + B * diff(b(x),x,2) + C * diff(c(x),x,2) == 0;
eqn2 = D * diff(d(x),x,2) - B * diff(b(x),x,2) - 2 * C * diff(c(x),x,2) - E * diff(e(x),x,2) == 0;
eqn3 = FF == (d(x) * b(x)) / a(x);
eqn4 = G == (d(x) * c(x)) / b(x);
eqn5 = H == d(x) * e(x);
eqns = [eqn1 eqn2 eqn3 eqn4 eqn5];
vars = [a(x); b(x); c(x); d(x); e(x)];
origVars = length(vars);
M = incidenceMatrix(eqns, vars);
[eqns, vars] = reduceDifferentialOrder(eqns, vars);
isLowIndexDAE(eqns,vars);
f = daeFunction(eqns,vars, A, B, C, D, E, FF, G, H);
A = 2.89e-9;   
B = 2.35e-9;  
C = 1.69e-9;  
D = 1.6e-8;    
E = 9.25e-9; 
FF = 6.24;      
G = 5.68e-5;  
H = 5.3e-8; 
F = @(t, Y, YP) f(t, Y, YP, A, B, C, D, E, FF, G, H);
vars;
y0est = [1.2650e4; 2.2748e4; 0.3724; 3.47; 1.5274e-8];
yp0est = zeros(5,1);
opt = odeset('RelTol', 10^(-3), 'AbsTol' , 10^(-3));
[y0, yp0] = decic(F, 0, y0est, [], yp0est, [], opt);
[tSol,ySol] = ode15i(F, [0, 7.21e-05], y0, yp0, opt);
for k = 1:origVars
  S{k} = char(vars(k));
end
plot(xSol,ySol(:,1),'r:',xSol,ySol(:,2),'k-.',xSol,ySol(:,3),'b--', xSol,ySol(:,4),'c-',xSol,ySol(:,5),'m-')
ylabel('Conc (mol/m^{3}')
xlabel ('Time (sec)')
legend('SO_2','HSO_{3}^{-}', 'SO_{3}^{2-}', 'H^+','OH^-','location','Best')

However I get the error message:

Index in position 1 exceeds array bounds (must not exceed 5).

Error in
symengine>@(x,in2,in3,param1,param2,param3,param4,param5,param6,param7,param8)[in3(6,:).*param1+in3(7,:).*param2+in3(8,:).*param3;-in3(7,:).*param2-in3(8,:).*param3.*2.0+in3(9,:).*param4-in3(10,:).*param5;param6-(in2(2,:).*in2(4,:))./in2(1,:);param7-(in2(3,:).*in2(4,:))./in2(2,:);param8-in2(4,:).*in2(5,:);in2(6,:)-in3(1,:);in2(7,:)-in3(2,:);in2(8,:)-in3(3,:);in2(9,:)-in3(4,:);in2(10,:)-in3(5,:)]

Error in WaterProfiles>@(t,Y,YP)f(t,Y,YP,A,B,C,D,E,FF,G,H)

Error in decic (line 66)
res = feval(odefun,t0,y0,yp0,varargin{:});

Error in WaterProfiles (line 45)
[y0, yp0] = decic(F, 0, y0est, [], yp0est, [], opt);

Debugging points me to line 45:

[y0, yp0] = decic(F, 0, y0est, [], yp0est, [], opt);

And when I check y0est, to me it does not exceed 5:

y0est = [1.2650e4; 2.2748e4; 0.3724; 3.47; 1.5274e-8];

What does "Index in position 1 exceeds array bounds (must not exceed 5)" mean?

[dpb]

y0est(1) --> 1.2650e4 >> 5

I don't have Symbolic TB so can't try to run the code to deicpher exactly what's going on but if it really is y0est being used as an indirect indexing array, the values therein as floating point aren't suitable.

What is decic?

[Dursman Mchabe]

Thanks a lot for the comment. I have also tried to use:

y0est = zeros(5,1);

And I still get the same error.

decic compute consistent initial conditions for ode15i. Its details are on:

https://www.mathworks.com/help/releases/R2018a/matlab/ref/decic.html

Answer 1

[Stephan]

Hi,

after you do:

[eqns, vars] = reduceDifferentialOrder(eqns, vars);

eqns and vars have a size of 10x1 sym both. So try the following:

y0est = [1.2650e4; 2.2748e4; 0.3724; 3.47; 1.5274e-8; 0; 0; 0; 0; 0];
yp0est = zeros(10,1);

It seems to work so far then - but then there a new Problem is coming up:

Undefined function or variable 'xSol'.

Error in alpha_z (line 32)
plot(xSol,ySol(:,1),'r:',xSol,ySol(:,2),'k-.',xSol,ySol(:,3),'b--', xSol,ySol(:,4),'c-',xSol,ySol(:,5),'m-')

Nowhere in your code xsol is defined - but in your function call of ode15i you define tsol - so replacing xsol by tsol gives a working code:

plot(tSol,ySol(:,1),'r:',tSol,ySol(:,2),'k-.',tSol,ySol(:,3),'b--', tSol,ySol(:,4),'c-',tSol,ySol(:,5),'m-')

I assume the problem you had is solved - even if the diagram does not look like something meaningful has been calculated - but I'm not a chemist ...

Best regards

Stephan

[Dursman Mchabe]

Thanks a lot Stephan. Indeed, my problem is solved. Thank you.