Array indices must be positive integers or logical values.

I created an ode45 function and wrote it like below

function rate = esb(t,C)

Q=[0.991 0.428 0.042 0.022 0.001 0.002 0.061 0.137 0.198]; % flowrate of the lake [m^3/s]

V=613411.183; % volume of the lake [m^3]

H=3.1; %depth of the lake [m]

u=[0.001106 0.000478 0.000047 0.000025 0.000001 0.000002 0.000068 0.000153 0.000221]; %[m/s]

k_a=3.93*(u.^0.5)/H^1.5; %reaeration rate constant calculated by O'Connor-Dobbins [1/day]

T=[15.3 25.1 24.0 28.6 26.2 23.4 16.0 14.3 8.5]; %temperature at [degree celcius]

k_s_20=0.15/H; %algal settling rate constant [m/day] at 20 celcius

theta_s=1.024;

k_s=k_s_20.*theta_s.^(T-20);

k_oa_20=0.21; %rate constant for the hydrolysis of organic N to Ammonium [1/day] at 20 celcius

theta_oa=1.047;

k_oa=k_oa_20.*theta_oa.^(T-20);

k_ai_20=0.5; %rate constant for the biological oxidation of NH4 to NO2 [1/day] at 20 celcius

theta_ai=1.083;

k_ai=k_ai_20.*theta_ai.^(T-20);

k_in_20=0.9; %rate constant for the biological oxidation of NO2 to NO3 [1/day] at 20 celcius

theta_in=1.047;

k_in=k_in_20.*theta_in.^(T-20);

k_1_20=0.6; %saturated algal growth rate constant [1/day] at 20 celcius

theta_1=1.068;

k_1=k_1_20.*theta_1.^(T-20);

k_2_20=0.1; %endegeneous respiration rate constant [1/day] at 20 celcius

theta_2=1.047;

k_2=k_2_20.*theta_2.^(T-20);

k_l_20=0.23; %carbonaceous deoxgeneration rate constant [1/day] at 20 celcius

theta_l=1047;

k_l=k_l_20.*theta_l.^(T-20);

F=0.5; %algal preference factor for ammonium

r_o1=1.104;

r_o2=0.955;

r_NA=0.063; %fraction of algal biomass that is nitrogen [mg Nitrogen/mg Algea]

r_oa=3.43; %fraction of algal biomass that is nitrogen [mg Nitrogen/mg Algea]

r_oi=1.14;

A_in=0.213; %inflow algea concentration [mg/L]

L_in=48.16; %inflow BOD concentration [mg/L]

N_org_in=0.3; %inflow organik nitrogen concentration [mg/L]

N_1_in=0.242; %inflow ammonium concentration [mg/L]

N_2_in=0.012; %inflow nitrite concentration [mg/L]

N_3_in=0.098; %inflow nitrate concentration [mg/L]

DO_in=8.2; %inflow dissolved oxygen concentration [mg/L]

DO_s=[10 8.2 8.4 7.7 8.1 8.5 9.9 10.2 11.7]; %saturated dissolved oxygen concentration [mg/L]

rate(1)=(Q(t)/V)*A_in-(Q(t)/V).*C(1)+k_1(t).*C(1)-k_2(t).*C(1)-k_s(t)*C(1);

rate(2)=(Q(t)/V)*L_in-(Q(t)/V)*C(2)-k_l(t)*C(2);

rate(3)=(Q(t)/V)*N_org_in-(Q(t)/V)*C(3)-k_oa(t)*C(3)+k_2(t)*r_NA*C(1);

rate(4)=(Q(t)/V)*N_1_in-(Q(t)/V)*C(4)+k_oa(t)*C(3)-k_ai(t)*C(4)-k_1(t)*F*r_NA*C(1);

rate(5)=(Q(t)/V)*N_2_in-(Q(t)/V)*C(5)+k_ai(t)*C(4)-k_in(t)*C(5);

rate(6)=(Q(t)/V)*N_3_in-(Q(t)/V)*C(6)+k_in(t)*C(5)-k_1(t)*(1-F)*r_NA*C(6);

rate(7)=(Q(t)/V)*DO_in-(Q(t)/V)*C(7)+k_a(t)*(DO_s(t)-C(7))+(k_1(t)*r_o1)*C(1)-(k_2(t)*r_o2)*C(1)-k_l(t)*C(2)-(k_ai(t)*r_oa)*C(4)-(k_in(t)*r_oi)*C(5);

rate=rate(:);

end

when I call this function in the script

ode45(@esb, [1 9], [7.903 56.7 3.533 0.533 0.062 0.072 8.4]);

legend('A','L','N_org','N_1','N_2','N_3','DO')

it doesn't work and give 'Array indices must be positive integers or logical values.' error. How can I solve this?

Answer 1

[J. Alex Lee]

In your function esb(), you are using the variable "t" as an index, e.g.,

Q(t)

but "t" is defined (via your argument structure) as the domain variable on which to solve the ODE.

This explains your error, but I don't think we have enough understanding of the underlying math problem to fix your code.

One possibility is that your parameters in your system of ODEs are all functions of the domain variable t, and you have hard-coded the values at some known values of "t" corresponding to integer values into your esb() function.

In this case, you will probably have to interpolate all of those parameters based on arbitrary values of "t" because your ODE solver chooses its own values. With minimal changes to your code, you'll have to change all of your X(t) to something like

t0 = 1:9;
Q_of_t = interp1(t0,Q,t);
% replace all instances of Q(t) with the above, same for k_a(t), etc.

If it were me, I would move all of these definitions out of the esb() function and supply them as function handles or something like that.