How to use lsqcurvefit to estimate parameters while using fsolve to solve an equation

Hi ebveryone,

On the attached code I am trying to use lsqcurvefit to estimate the following constants:

% Constants
KSO2    = Constant(1);
KHSO3   = Constant(2);
KCO2    = Constant(3);
KHCO3   = Constant(4);
Kw      = Constant(5);

whilst using fsolve to solve an equation. Some I get an error message:

Error using pHCalcFeb25>pHCalculation
Too many output arguments.

Error in lsqcurvefit (line 213)
            initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});

Error in pHCalcFeb25 (line 67)
    Constant=lsqcurvefit(@pHCalculation,Constant0,t,pH);

Caused by:
    Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.
 
>> 

How can I resolve this error. Please help

Answer 1

[Walter Roberson]

Your function calculates ph as a vector of 31 values. However, your function is declared as if it has no return values at all. At the very least you need to change your function declaration to return ph .

It was easier to debug by splitting your function into two files, which I attach.

[Dursman Mchabe]

Thanks a lot Walter. I made a mistake on line 172:

Cfit = kinetics(Constant, tv);

It was supposed to be:

Cfit = pHCalculation1(Constant, tv);

This solves the previous error, yet it now gives:

Not enough input arguments.

Error in pHCalculation1 (line 15)
    sol(i)    = fsolve(@(x)  x + 2.* Ca_total(i) - ((S_total(i).* KSO2.*x)/(x.^2 + KSO2.*x + KSO2*KHSO3))- 2.*((S_total(i).*KSO2*KHSO3)/(x.^2 ...

Error in pHCalc1Feb25 (line 172)
Cfit = pHCalculation1(Constant, tv);

[Dursman Mchabe]

When I run the code with just dummy values, I get success e.g:

function pHCalculation()
% Constant Parameters
KSO2    = 6.24;
KHSO3   = 5.68e-5;
KCO2    = 1.7e-3;
KHCO3   = 6.55e-8;
Kw      = 5.3e-8;
CH0     = 7.41e-06;

t       = [0
        6000
        12000
        18000
        24000
        30000
        36000
        42000
        48000
        54000
        60000
        66000
        72000
        78000
        84000
        90000
        96000
        102000
        108000
        114000
        120000
        126000
        132000
        138000
        144000
        150000
        156000
        162000
        168000
        174000
        180000];

Ca_total   = [0.04879
6.650241034
9.984349357
12.80945024
15.8072316
18.85920165
21.81901476
24.15458161
23.99797436
22.9179308
22.08000295
21.3980477
20.84506067
20.40448911
20.06253803
19.80636092
19.623386
19.50098984
19.42639218
19.3868172
19.37002007
19.36526318
19.36471126
19.37750389
19.43991156
19.51435862
19.59106428
19.66791631
19.74434773
19.82019925
19.89542823
];

S_total   = [0
3.228125432
3.039568727
2.353698976
1.924485852
1.669257876
1.577120575
1.848544338
3.227063
5.416789376
7.795831728
10.27649137
12.83100297
15.43762835
18.07548351
20.72364934
23.36086852
25.96529992
28.51423028
30.98359907
33.34665967
35.5684028
37.57358167
39.0779715
39.84164253
40.20449967
40.43153695
40.61599762
40.78655902
40.95174085
41.11405033];

C_total = [0.04879
        6.89871101
        13.90837796
        20.90460441
        27.81415915
        34.59918462
        41.11877195
        46.61855348
        48.32509945
        47.89047761
        47.18901391
        46.2708157
        45.16329308
        43.89007529
        42.47196783
        40.92750491
        39.27332492
        37.52445179
        35.69453309
        33.79609251
        31.84091599
        29.84097095
        27.8120307
        25.79309152
        23.85163
        22.01053317
        20.26841487
        18.62064555
        17.06221895
        15.58824986
        14.19408257];
    

pH      = [8.13
        7.76
        7.43
        7.45
        7.45
        7.00
        4.39
        3.61
        3.21
        2.88
        2.88
        2.64
        2.63
        2.62
        2.74
        2.86
        2.86
        2.68
        2.74
        2.8
        2.64
        2.79
        2.79
        2.95
        3.05
        3.07
        3.12
        3.06
        3.06
        2.95
        3.00];
    
for i = 1:length(t)
      sol(i)    = fsolve(@(x)  x + 2.* Ca_total(i) - ((S_total(i).* KSO2.*x)/(x.^2 + KSO2.*x + KSO2*KHSO3))- 2.*((S_total(i).*KSO2*KHSO3)/(x.^2 ...
              + KSO2.*x + KSO2*KHSO3))-((C_total(i).*KCO2.*x)/(x.^2 + KCO2.*x + KCO2*KHCO3))- 2.*((C_total(i).*KCO2.*KHCO3)/(x.^2 ...
              + KCO2.*x + KCO2*KHCO3))- Kw./x, CH0);
end
ph = 3 - log10(sol);
plot(t,ph,'r--')
hold on
plot (t,pH,'kd')
xlabel('Time (sec)')
ylabel('pH')
hold off
legend('pH-Mod','pH-Exp')
end

[Walter Roberson]

What is phCalculation1 ?

[Walter Roberson]

Your phCalculation routine does not use the value of t at any point, only the length of t. If length(t) is greater than length() of Ca_total and the other arrays, then you would get an error attempting to access those arrays; otherwise you would get the first length(t) entries returned instead of the full 31. linspace() with no size parameter uses 100 by default, which is more than the number of entries in the array.

It looks to me as if what you would want is to interpolate the results at those times. But on the other hand you plot() with t (length 31) as the axes rather than with tv (length 100).

Perhaps

Cfit = pHCalculation(Constant, t, Ca_total, S_total, C_total);

[Dursman Mchabe]

So sorry, I already had the files @pHCalculation and @pHCalcFeb25 in the folder, and when you sent those two m-files saved them as pHCalculation1 and pHCalc1Feb25, just to quickly test them.

Apologies, I have now moved them to a different folder, and saved them as pHCalculation and pHCalcFeb25.

Sorry for causing confusion.

[Dursman Mchabe]

I Actually need a tv of length 31. I will try:

Cfit = pHCalculation(Constant, t, Ca_total, S_total, C_total);

or

tv = linspace(min(t), max(t),31);

[Dursman Mchabe]

Thank you , thank you , thank you. It works perfectly.