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

Dursman Mchabe 님이 질문을 제출함. 24 Feb 2019

최근 활동 Dursman Mchabe 님이 댓글을 추가함. 25 Feb 2019

Walter Roberson 님이 답변을 채택함.

조회 수: 3(최근 30일)

pHCalcFeb25.m

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

댓글 수: 0

로그인 to comment.

답변 수: 1

링크

Walter Roberson 님의 답변 25 Feb 2019

채택된 답변

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.

댓글 수: 7

표시 이전 댓글 수: 4

Dursman Mchabe 25 Feb 2019

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

Dursman Mchabe 25 Feb 2019

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 25 Feb 2019

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

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

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

Dursman Mchabe 25 Feb 2019

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

로그인 to comment.

이 질문에 답변하려면 로그인을(를) 수행하십시오.

이 질문을 팔로우합니다.

Dursman Mchabe (view profile)