Estimation of parameters and initial conditions using lsqcurvefit gives an error message :"Index Exceeds Array Bounds"

1 view (last 30 days)
Dursman Mchabe
Dursman Mchabe on 30 Sep 2018
Commented: Dursman Mchabe on 1 Oct 2018
Hi everyone, I tried different approaches to emulate Dave Black's method:
https://www.mathworks.com/matlabcentral/answers/329901-solving-coupled-differential-equations
But now, I get the error message:
Index exceeds array bounds.
Error in sym/subsref (line 859)
R_tilde = builtin('subsref',L_tilde,Idx);
Error in EstimationOfParametersAndInitialConditionsUsingLSQCURVEFIT1/ode15ifun (line 75)
eqn1 = diff(b(5),t)== F_rate/ V_Headspace * ( CSO2_in - b(5)) - (b(5) * R * T - HSO2 * ((b(7) * b(12)^2)/(b(12)^2 + KSO2
* b(12) + KSO2 * KHSO3))) / ((1/b(1)) + 1/(b(2) * (((DSO2 * (HSO2 * b(5) * R * T - ((b(7) * b(12)^2)/(b(12)^2 + KSO2 *
b(12) + KSO2 * KHSO3))) + DHSO3 * (((KSO2 * HSO2 * b(5) * R * T / b(12)) - ((b(7) * KSO2 * b(12))/(b(12)^2 + KSO2 * b(12)
+ KSO2 * KHSO3)))) + DSO32 * ((((KSO2 * KHSO3 * HSO2 * b(5) * R * T/b(12)^2) - ((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 *
b(12) + KSO2 * KHSO3)))))) / (DSO2 * (HSO2 * b(5) * R * T - ((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 *
KHSO3))))))));
Error in lsqcurvefit (line 213)
initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});
Error in EstimationOfParametersAndInitialConditionsUsingLSQCURVEFIT1 (line 155)
B = lsqcurvefit(@ode15ifun, eqns0, ExperimentalTime, ExperimentalConcentrations, zeros(size(eqns0)), inf(size(eqns0)))
Caused by:
Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.
My current code is:
function EstimationOfParametersAndInitialConditionsUsingLSQCURVEFIT1
ExperimentalTime = [600
1200
1800
2400
3000
10200
17400
24600
31800
39000
46200
53400
60600
67800
75000
82200
89400
96600
103800
111000
118200
125400
132600
139800
147000
154200
161400
168600
175800
183000];
ExperimentalConcentrations = [ 0.00E+00 6.97E-08 3.34E-01 4.88E-06 4.88E-06 4.86E+01 0.00E+00 7.41E-06
1.72E-02 8.76E-08 6.79E-01 4.88E-06 4.88E-06 4.92E+01 0.00E+00 9.33E-06
1.71E-02 1.18E-07 1.04E+00 5.11E-06 5.11E-06 4.94E+01 0.00E+00 1.20E-05
1.70E-02 1.03E-07 1.41E+00 5.11E-06 5.11E-06 4.96E+01 0.00E+00 1.05E-05
1.68E-02 1.77E-07 6.00E+00 5.35E-06 5.35E-06 3.78E+01 6.30E+00 1.74E-05
1.74E-02 3.65E-07 1.07E+01 5.35E-06 5.35E-06 3.09E+01 9.33E+00 3.72E-05
1.65E-02 3.50E-07 1.56E+01 5.36E-06 5.36E-06 2.67E+01 1.13E+01 3.55E-05
1.67E-02 8.83E-07 2.07E+01 5.36E-06 5.36E-06 1.98E+01 1.43E+01 1.00E-04
1.67E-02 5.01E-06 2.59E+01 5.07E-06 5.07E-06 1.14E+01 1.79E+01 4.07E-02
1.70E-02 5.06E-06 3.14E+01 5.07E-06 5.07E-06 9.58E+00 1.91E+01 2.45E-01
1.80E-02 4.93E-06 3.67E+01 4.93E-06 4.93E-06 5.65E+00 2.09E+01 6.17E-01
2.15E-02 5.29E-06 4.18E+01 5.30E-06 5.30E-06 3.32E+00 2.21E+01 1.32E+00
2.68E-02 5.30E-06 4.66E+01 5.30E-06 5.30E-06 0.00E+00 2.37E+01 2.29E+00
3.07E-02 5.30E-06 5.09E+01 5.30E-06 5.30E-06 0.00E+00 2.40E+01 2.34E+00
3.33E-02 9.45E-06 5.51E+01 9.46E-06 9.46E-06 0.00E+00 2.43E+01 2.40E+00
3.68E-02 9.45E-06 5.90E+01 9.46E-06 9.46E-06 0.00E+00 2.47E+01 1.82E+00
4.14E-02 1.13E-05 6.23E+01 1.13E-05 1.13E-05 0.00E+00 2.50E+01 1.38E+00
4.48E-02 1.15E-05 6.56E+01 1.15E-05 1.15E-05 0.00E+00 2.54E+01 2.09E+00
4.87E-02 1.10E-05 6.87E+01 1.10E-05 1.10E-05 0.00E+00 2.58E+01 1.82E+00
5.19E-02 1.10E-05 7.12E+01 1.10E-05 1.10E-05 0.00E+00 2.62E+01 1.58E+00
5.47E-02 1.04E-05 7.36E+01 1.04E-05 1.04E-05 0.00E+00 2.67E+01 2.29E+00
5.75E-02 1.04E-05 7.59E+01 1.04E-05 1.04E-05 0.00E+00 2.71E+01 1.62E+00
5.97E-02 1.04E-05 7.78E+01 1.04E-05 1.04E-05 0.00E+00 2.76E+01 1.12E+00
6.17E-02 1.03E-05 7.95E+01 1.03E-05 1.03E-05 0.00E+00 2.81E+01 8.91E-01
6.33E-02 1.04E-05 8.11E+01 1.04E-05 1.04E-05 0.00E+00 2.86E+01 8.51E-01
6.44E-02 1.05E-05 8.24E+01 1.06E-05 1.06E-05 0.00E+00 2.90E+01 7.59E-01
6.48E-02 1.13E-05 8.24E+01 1.13E-05 1.13E-05 0.00E+00 2.90E+01 8.71E-01
6.53E-02 1.13E-05 8.23E+01 1.13E-05 1.13E-05 0.00E+00 2.90E+01 1.12E+00
6.57E-02 1.13E-05 8.21E+01 1.13E-05 1.13E-05 0.00E+00 2.90E+01 1.00E+00
6.59E-02 1.03E-05 8.20E+01 1.03E-05 1.03E-05 0.00E+00 2.90E+01 8.51E-01];
%-------------------------------------------------------------------------------------------------------------------------------------
function eqns = ode15ifun(b,t,V_Headspace, F_rate, CSO2_in, CCO2_in, R, T, HSO2, DCa2, DSO2, DHSO3, DSO32, kLa_CO2, HCO2, DCO2, DHCO3, DCO32, KSPCaSO3, BETCaSO3, BETCaCO3, MWCaCO3, KCO2, KHCO3, KSO2, KHSO3, Kw, CCa2)
% b(1:4) = Parameters, b(5:12) = Initial Conditions
% MAPPING: b(1) = kga, b(2) = kLa_SO2, b(3) = ktot, b(4) = Kad,
% b(5)= CSO2_gas(t), b(6)= CCO2_gas(t), b(7)= S_total(t), b(8)= C_total(t), b(9)= Ca2_total(t),b(10)= CCaCO3(t),b(11)= CCaSO3(t),b(12)= CH(t)
syms b V_Headspace F_rate CSO2_in CCO2_in R T HSO2 DCa2 DSO2 DHSO3 DSO32 kLa_CO2 HCO2 DCO2 DHCO3 DCO32 KSPCaSO3 BETCaSO3 BETCaCO3 MWCaCO3 KCO2 KHCO3 KSO2 KHSO3 Kw CCa2
eqn1 = diff(b(5),t)== F_rate/ V_Headspace * ( CSO2_in - b(5)) - (b(5) * R * T - HSO2 * ((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))) / ((1/b(1)) + 1/(b(2) * (((DSO2 * (HSO2 * b(5) * R * T - ((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))) + DHSO3 * (((KSO2 * HSO2 * b(5) * R * T / b(12)) - ((b(7) * KSO2 * b(12))/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3)))) + DSO32 * ((((KSO2 * KHSO3 * HSO2 * b(5) * R * T/b(12)^2) - ((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3)))))) / (DSO2 * (HSO2 * b(5) * R * T - ((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))))))));
eqn2 = diff(b(6) ,t)== F_rate/ V_Headspace * ( CCO2_in - b(6)) - (kLa_CO2 * ((((DCO2*(HCO2 * b(6) * R * T - ((b(8) * b(12)^2)/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3))) + DHCO3*(((KCO2 * HCO2 * b(6) * R * T/b(12)) - ((C_total (t) * KCO2 * b(12))/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3)))) + DCO32 * ((((KCO2 * KHCO3 * HCO2 * b(6) * R * T/b(12)^2) - ((b(8) * KCO2 * KHCO3)/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3)))))) / (DCO2 * (HCO2 * b(6) * R * T - ((C_total (t) * b(12)^2)/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3)))))))*(b(6) * R * T/HCO2 - ((b(8) * b(12)^2)/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3))));
eqn3 = diff(b(7),t)== (b(5) * R * T - HSO2*((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))) / ((1/b(1)) + 1/(b(2) * (((DSO2*(HSO2* b(5) * R * T - ((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))) + DHSO3*(((KSO2 * HSO2 * b(5) * R * T/b(12)) - ((b(7) * KSO2 * b(12))/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3)))) + DSO32 * ((((KSO2 * KHSO3 *HSO2* b(5) * R * T/b(12)^2) - ((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3)))))) / (DSO2 * (HSO2 * b(5) * R * T - ((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))))))) - 0.162 * exp(-5153/T) * (((CCa2 * ((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))) / KSPCaSO3) - 1)^2 * (BETCaSO3/((CCa2 * (((((KSO2 * KHSO3 *HSO2* b(5) * R * T/b(12)^2) - ((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3)))))) / (DSO2 * (HSO2 * b(5) * R * T - ((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))))))))/ KSPCaSO3;
eqn4 = diff(b(8),t)== (kLa_CO2 * ((((DCO2*(HCO2 * b(6) * R * T - ((b(8) * b(12)^2)/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3))) + DHCO3*(((KCO2 * HCO2 * b(6) * R * T/b(12)) - ((C_total (t) * KCO2 * b(12))/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3)))) + DCO32 * ((((KCO2 * KHCO3 * HCO2 * b(6) * R * T/b(12)^2) - ((b(8) * KCO2 * KHCO3)/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3)))))) / (DCO2 * (HCO2 * b(6) * R * T - ((C_total (t) * b(12)^2)/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3)))))))*(b(6) * R * T /HCO2 - ((b(8) * b(12)^2)/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3)))) + (b(3) * BETCaCO3 * MWCaCO3 * b(10) * b(12) *(1 - (b(4) * b(12))/(1 + b(4) * b(12))));
eqn5 = diff(b(9),t)== (-1) * (b(3) * BETCaCO3 * MWCaCO3 * b(10) * b(12) *(1 - (b(4) * b(12))/(1 + b(4) * b(12)))) - (0.162 * exp(-5153/T) * (((CCa2 * ((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))) / KSPCaSO3) - 1)^2 * (BETCaSO3/((CCa2 * (((((KSO2 * KHSO3 *HSO2* b(5) * R * T/b(12)^2) - ((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3)))))) / (DSO2 * (HSO2 * b(5) * R * T - ((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))))))))/ KSPCaSO3;
eqn6 = diff(b(10),t)== (-1) * (b(3) * BETCaCO3 * MWCaCO3 * b(10) * b(12) * (1 - (b(4) * b(12))/(1 + b(4) * b(12))));
eqn7 = diff(b(11),t)== 0.162 * exp(-5153/T) * ((((CCa2 * ((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))) / KSPCaSO3) - 1)^2 * (BETCaSO3/((CCa2 * (((((KSO2 * KHSO3 *HSO2* b(5) * R * T/b(12)^2) - ((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3)))))) / (DSO2 * (HSO2 * b(5) * R * T - ((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))))))))/ KSPCaSO3;
eqn8 = diff(b(12),t)== b(12) + 2 * CCa2 - ((b(7) *KSO2 * b(12))/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))- 2*((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3)) - ((b(8) * KCO2 * b(12))/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3)) - 2 * ((b(8) * KCO2 * KHCO3)/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3))- Kw / b(12);
eqns = [eqn1 eqn2 eqn3 eqn4 eqn5 eqn6 eqn7 eqn8];
vars = [b(5); b(6); b(7); b(8); b(9); b(10); b(11); b(12)];
origVars = length(vars);
M = incidenceMatrix(eqns, vars);
[eqns, vars] = reduceDifferentialOrder(eqns, vars);
isLowIndexDAE(eqns,vars);
f = daeFunction(eqns,vars, V_Headspace, F_rate, CSO2_in, CCO2_in, CCa2, R, T, HSO2, b(1), DCa2, DSO2, DHSO3, DSO32, b(2), kLa_CO2, HCO2, DCO2, DHCO3, DCO32, KSPCaSO3, BETCaSO3, b(3), BETCaCO3, MWCaCO3, b(4), KCO2, KHCO3, KSO2, KHSO3, Kw);
V_Headspace = 1.5e-6;
F_rate = 1.66667e-5;
CSO2_in = 6.51332e-2;
CCO2_in = 0;
CCa2 = 10;
R = 8.314;
T = 323.15;
HSO2 = 149;
b(1) = 4.14e-6;
DCa2 = 1.39e-9;
DSO2 = 2.89e-9;
DHSO3 = 2.89e-9;
DSO32 = 2.89e-9;
b(2) = 8.4e-4;
kLa_CO2 = 9.598e-4;
HCO2 = 5.15e3;
DCO2 = 3.53e-9;
DHCO3 = 2.89e-9;
DCO32 = 2.89e-9;
KSPCaSO3 = 1.07e-7;
BETCaSO3 = 10;
b(3) = 8.825e-3;
BETCaCO3 = 12.54;
MWCaCO3 = 100.0869;
b(4) = 0.84;
KCO2 = 1.7e-3;
KHCO3 = 6.55e-8;
KSO2 = 6.24;
KHSO3 = 5.68e-5;
Kw = 5.3e-8;
F = @(t, Y, YP) f(t, Y, YP, V_Headspace, F_rate, CSO2_in, CCO2_in, CCa2, R, T, HSO2, b(1), DCa2, DSO2, DHSO3, DSO32, b(2), kLa_CO2, HCO2, DCO2, DHSO3, DSO32, KSPCaSO3, BETCaSO3, b(3), BETCaCO3, MWCaCO3, b(4), KCO2, KHCO3, KSO2, KHSO3, Kw);
vars;
%b0est = [0; 0; 3.3351e-1; 4.879e-6; 4.879e-6; 4.862e+1; 0; 7.413e-6];
b0est = [1;1;1;1;1;1;1;1;];
bp0est = zeros(8,1);
opt = odeset('RelTol', 10.0^(-7), 'AbsTol' , 10.0^(-7));
[b0, bp0] = decic(F, 0, b(5:12), [], bp0est, [], opt);
[tSol,eqns] = ode15i(F, [600, 183000], b0, bp0, opt);
for k = 1:origVars
S{k} = char(vars(k));
end
end
eqns0 = [4.14e-6; 8.4e-4; 8.825e-3; 0.84; 0; 0; 3.3351e-1; 4.879e-6; 4.879e-6; 4.862e+1; 0; 7.413e-6];
B = lsqcurvefit(@ode15ifun, eqns0, ExperimentalTime, ExperimentalConcentrations, zeros(size(eqns0)), inf(size(eqns0)))
fprintf(1, '\n\tParameters:\n\t\tkga = %11.3E\n\t\tkLa = %11.3E\n\t\tktot = %11.3E\n\t\tKad = %11.3E\n\tInitial Conditions:\n\t\tCSO2gas = %11.3E\n\t\tCCO2gas = %11.3E\n\t\tStotal = %11.3E\n\t\tCtotal = %11.3E\n\t\tCatotal = %11.3E\n\t\tCCaCO3 = %11.3E\n\t\tCCaSO3 = %11.3E\n\t\tCH = %11.3E\n', B);
F = @(t, Y, YP) f(t, B, BP, V_Headspace, F_rate, CSO2_in, CCO2_in, CCa2, R, T, HSO2, B(1), DCa2, DSO2, DHSO3, DSO32, B(2), kLa_CO2, HCO2, DCO2, DHSO3, DSO32, KSPCaSO3, BETCaSO3, B(3), BETCaCO3, MWCaCO3, B(4), KCO2, KHCO3, KSO2, KHSO3, Kw);
[tSol,eqns] = ode15i(F, [600, 183000], B(5:12), bp0, opt);
figure(1)
plot(ExperimentalTime, ExperimentalConcentrations, 'p')
hold on
plot(tSol, eqns, '--')
hold off
grid
legend('Exp-SO_{2-gas}','Exp-CO_{2-gas}','Exp-S_{total}','Exp-C_{total}','Exp-Ca^{2+}_{total}','Exp-CaCO_{3}','Exp-CaSO_{3}.^{1}/_{2}H_{2}O','Exp-H^{+}', 'Mod-SO_{2-gas}','Mod-CO_{2-gas}','Mod-S_{total}','Mod-C_{total}','Mod-Ca^{2+}_{total}','Mod-CaCO_{3}','Mod-CaSO_{3}.^{1}/_{2}H_{2}O','Mod-H^{+}', 'Location','best')
end
What can I try?

Sign in to answer this question.

Answers (1)

Torsten
Torsten on 1 Oct 2018

In

[tSol,eqns] = ode15i(F, [600, 183000], b0, bp0, opt);

you must specify the complete vector "ExperimentalTime", not only [600, 183000].

Best wishes

Torsten.

  1 Comment
Dursman Mchabe
Dursman Mchabe on 1 Oct 2018
Thanks a lot for the comment. I have specified the complete vector "ExperimentalTime", however the error does not go away. I am currently inspecting the code, and thinking of the next move.

Sign in to comment.

Categories

Find more on Equation Solving in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!