Index out of bound in ODE Solver

So I have a basic combustion model where the dependent variables are highly dependent on each other:

function ydot = ss_combustion_counter(t,y)

ydot(1) = theta;

ydot(2) = -(1/lambda_solid)*(r_char_oxi*M_C*h_char_oxi + r_p1*h_p1) + (1/lambda_solid)*h_sg*Apvp*(Ts - Tf)... + (1/lambda_solid)*4*167*(295-Ts)/0.125;

ydot(3) = h_sg*Apvp*(Ts - Tf) + 4*1.256*(295-Tf)/0.125 + (- r_tar*M_tar*h_tar - r_ch4*M_CH4*h_ch4 - r_co*M_CO*h_co - r_h2*M_H2*h_h2);

ydot(4) = (r_p1*h2o_p1 + (r_tar*0.761*M_H2O) + (r_ch4*2*M_H2O) + (r_h2*M_H2O) + (r_char_oxi*(0.2526/2)*M_H2O) - y(4)*ydot(13))/y(13);

ydot(5) = (r_p1*co_p1 + (r_tar*M_CO) + (r_ch4*M_CO) - (r_co*M_CO) + (r_char_oxi*(2-2*0.8013-0.0237+0.2526/2)*M_CO) - y(5)*ydot(13))/y(13);

ydot(6) = (r_p1*co2_p1 + (r_co*M_CO2) + (r_char_oxi*(2*0.8013+0.0237-0.2426/2+1)*M_CO2) - y(6)*ydot(13))/y(13);

ydot(7) = (r_p1*h2_p1 - (r_h2*M_H2) - y(7)*ydot(13))/y(13);

ydot(8) = (r_p1*ch4_p1 - (r_ch4*M_CH4) - y(8)*ydot(13))/y(13);

ydot(9) = (r_p1*tar_p1 - r_tar*M_tar - y(9)*ydot(9))/y(13);

ydot(10) = ((-r_tar*0.867*M_O2) - (r_ch4*1.5*M_O2) - (r_co*0.5*M_O2) - (r_h2*0.5*M_O2) - (r_char_oxi*0.8013*M_O2) - y(10)*ydot(13))/y(13);

ydot(11) = 0;

ydot(12) = ydot(4) + ydot(5) + ydot(6) + ydot(7) + ydot(8) + ydot(9) + ydot(10) + ydot(11);

ydot(13) = -(y(13)/y(12)) * ydot(12) + (1/y(12))*r_p1*h2o_p1 + (1/y(12))*gas_reaction_sum;

ydot = ydot';

Other physical parameters are defined in a separate main file. When I run it, I get an error saying:

Attempted to access ydot(13); index out of bounds because numel(ydot)=3

So I can see that it's because I'm using using ydot(13) term before it's defined in the ydot(4) expression. But I'm a little confused why this would be an issue since MATLAB allows you do use higher value such as y(13) in the ydot(4) expression without too much trouble. Since ydot(4)~ydot(11) and ydot(12), ydot(13) expressions are interdependent on each other, I'm not sure how to get around this issue.

If anyone has a good idea, it'd greatly help with my project. Thank You! - Charlie

Answer 1

[Star Strider]

I didn’t run your code, but since I couldn’t find any element of ydot in a denominator, I suggest you preallocate ydot by putting:

ydot = zeros(13,1);

as the first line in ss_combustion_counter.

See also the section on ‘Passing Extra Parameters’ in the ODE solver documentation if you need that. It’s not obvious to me how they’re getting into ss_combustion_counter now, but I may be missing something.

[Charlie]

Thanks a lot. This let my program run without running into the index issue. One question I have is since I defined ydot = zeros(13,1);, wouldn't the program recognize ydot(13) = 0, not what I defined above which is:

ydot(13) = -(y(13)/y(12)) * ydot(12) + (1/y(12))*r_p1*h2o_p1 + (1/y(12))*gas_reaction_sum;

Thanks~

[Star Strider]

My pleasure!

Defining ydot = zeros(13,1) at the beginning of your function simply preallocates (creates space in memory for) the entire 13-element ydot vector by creating a (13x1) column vector of zeros. The values for ydot(1) through ydot(13) are then calculated and returned to the ODE solver as a 13-element column vector with the values calculated as you defined them.