ode45 function array error

1 view (last 30 days)
Carey n'eville
Carey n'eville on 13 Jan 2021
Edited: James Tursa on 13 Jan 2021
Hi friends there occur when I run my code. Could you help me please for solving my problem? Code and error are given below:
THE CODE:
global mu_H Y_H Ks X_BH teta mu_Hmax T Ss_in
mu_Hmax=0.15; %days^-
mu_Hmax=mu_Hmax*((1.072)^(T-20)); %day^-1
T=22.6; %centigrade
X_BH=2295; %mg/l
teta=1.25; %day
Y_H=0.6; %g/g
Ks=250; %mg/l
%initial condition is Ss_in=1127 mg/l
Ss_in=1127; %mg/l
% Ss_in range is taken as +- %50 of Ss_,in value given
Sinr=round(Ss_in*0.5/10)*10:10:round(Ss_in*1.5/10)*10;
tspan=0:140;
tm=[0 30 32 34 36 50 52 57 59 61 63 65 72 75 76 79 83 85 87 89 96 99 106 109 116 118 121 124 125 126 127 131 133 134 135 136 137 140];
Ss_in=ones(length(tspan),1)*Ss_in;
res=zeros(length(tm),length(Sinr));
for i=1:length(tm)
for j=1:length(Sinr)
% update Sin value at each step for each value in the Sin range
Ss_in(tm(i):end)=Sinr(j);
ode45(@substrate,tspan,Ss_in)
% calculate error the specific value of S at that time step
res(i,j)=sqrt((Ss(tm(i)+1)-Sm(i))^2);
end
% find the location of minimum error
[minres, lmin]=min(res(i,:));
% update Sin vector with the value that gives minimum error
% at each time step
Ss_in(tm(i):end)=Sinr(lmin);
end
xlabel('Time (day)')
ylabel('Concentrations (mg/L)')
legend ('Substrate')
clear i j
Sp=zeros(length(tm),1);
for i=1:length(tm)
Sp(i)=Ss(tm(i)+1);
end
RMSE=sqrt(sum((Sp-Sm).^2));
%substrate function is defined in order to calculate dSs/dt, predictions
function Derivative = substrate(t,Ss)
global mu_H Y_H Ks X_BH teta Ss_in
Derivative =-mu_H/Y_H*(Ss/(Ks+Ss))*X_BH+(1/teta)*(Ss_in(round(t+1))-Ss); %dSs/
Derivative=Derivative(:);
end
ERROR:
Array indices must be positive integers or logical values.
Error in editTPHW3P3S1 (line 23)
Ss_in(tm(i):end)=Sinr(j);

Answers (1)

James Tursa
James Tursa on 13 Jan 2021
Edited: James Tursa on 13 Jan 2021
You have tm(1) = 0, so using tm(1) as an index is invalid:
tm=[0 30 etc...
:
for i=1:length(tm)
:
Ss_in(tm(i):end)=Sinr(j);

Categories

Find more on Programming in Help Center and File Exchange

Products


Release

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!