ode45 function array error

2 views (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 Particle & Nuclear Physics 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!