Fitting model to titration data

25 次查看（过去 30 天）
Rafael Ibáñez
Rafael Ibáñez 2019-4-12
评论： Rafael Ibáñez 2019-4-13
I'm trying to fit a model of titration to experimental data. The function "calcpH" calculate a pH value for each experimental point and I'm trying to obtain the value of Ca (and eventually Ka) that fit best the model to the experimental values.
When I uncomment the lines of the fitting procedure I get the following error:
Error using ajuste_listado_1>sseval
Too many output arguments.
Error in ajuste_listado_1>@(x)sseval(x,vbr,pHr)
Error in fminsearch (line 200)
fv(:,1) = funfcn(x,varargin{:});
What is going wrong ????
THANK YOU
CODE OF THE SCRIPT
clear
clc
pHr=[3.15,3.30,3.48,3.65,3.73,3.86,3.93,4.03,4.10,4.20,4.28,4.35,4.44,4.51,4.59,4.66,4.75,4.84,4.94,5.03,5.18,5.34,5.54,5.91,6.10,6.40,8.66,9.32,9.88,10.16,10.39,10.49,10.64,10.97,11.16,11.28,11.37,11.45,11.51,11.56,11.60,11.67];
vbr=[0.00 ,0.20 ,0.40 ,0.60 ,0.80 ,1.00 ,1.20 ,1.40 ,1.60 ,1.80 ,2.00 ,2.20 ,2.40 ,2.60 ,2.80 ,3.00 ,3.20 ,3.40 ,3.60 ,3.80 ,4.00 ,4.20 ,4.40 ,4.60 ,4.70 ,4.80 ,4.90 ,5.00 ,5.10 ,5.20 ,5.30 ,5.40 ,5.50 ,6.00 ,6.50 ,7.00 ,7.50 ,8.00 ,8.50 ,9.00 ,9.50 ,10.00];
Ka=1.8e-5;
Ca=0.01
Va=50;
Cb=0.1;
Vb=0:0.2:10;
Kw=1e-14;
% FITTING
% fun=@(x)sseval(x,vbr,pHr);
% bestx = fminsearch(fun,x0)
% Ca = bestx(1);
x=zeros(1)
pH = calcpH(x,Ka,Ca,Kw,Cb,Va,vbr);
plot(vbr,pHr,'o',vbr,pH,'x')
title('Titration of acetic acid')
xlabel('Vol NaOH [mL]')
ylabel('pH')
function sseval(x,vbr,pHr)
Ca=x(1)
pH = calcpH(x,Ka,Ca,Kw,Cb,Va,vbr)
sse = sum((pHr - pH').^2)
end
function pH=calcpH(x,Ka,Ca,Kw,Cb,Va,vbr)
VT=Va+vbr;
for i=1:length(vbr)
syms xh
h=vpasolve(Cb*vbr(i)/VT(i)-(Ca*Va/VT(i))*Ka/(Ka+xh)-Kw/xh+xh==0,xh);
pH(i)= -log10(h(3));
end
end

请先登录，再回答此问题。

采纳的回答

Rafael Ibáñez
Rafael Ibáñez 2019-4-13
Thank you very much

更多回答（1 个）

David Wilson
David Wilson 2019-4-13
编辑：David Wilson 2019-4-13
For one thing, you don't assign an answer to function sseval. That's not going to help things.
Here's my solution:
I've re-written your internal function and removed the numerical root finder. It's only a cubic, so we should just use roots. HOWEVER you have assumed you only get one physcially sensible (i.e. positive) root and use that. Are you sure that's OK? (I followed your approach, so I'm counting on that.)
function pH=calcpH(Ca,Ka,Kw,Cb,Va,vbr)
VT=Va+vbr;
a = Cb*vbr./VT;
b = Ca*Va*Ka./VT;
c = Ka;
d = Kw;
pH = NaN(size(vbr)); % placeholder
for i=1:length(vbr)
coeff3 = [1, (a(i)+c), -(b(i)+d-a(i)*c), -c*d];
h = sort(roots(coeff3)); % assume last is the only (sensible) positive
h = h(end);
pH(i,1)= -log10(h);
end
end
Now I wrap the objective function around that. I'm using the 2-norm, nice and well-conditioned. The output of this function is ideally small, or better still zero.
function sse = sseval(Ca,vbr,pHr, ...
Ka,Kw,Cb,Va)
pH=calcpH(Ca,Ka,Kw,Cb,Va,vbr);
sse = norm(pHr-pH); % sum((pHr - pH).^2);
end
Note that you need to pass through all your constants etc correctly.
Now we are ready to call the optimiser from the top script. I like to make sure everything is in columns (your original version had a mixture which was just asking for trouble.) You need a guess for Ca.
pHr = pHr(:); vbr = vbr(:); % ensure columns
fun=@(x) sseval(x, vbr,pHr, Ka,Kw,Cb,Va);
Ca_guess = 0.01;
optns = optimset('fminsearch');
[Ca, fval, exitflag] = fminsearch(fun,Ca_guess, optns)
%Ca = fminunc(fun,Ca_guess)
pH = calcpH(Ca,Ka,Kw,Cb,Va,vbr);
plot(vbr,pHr,'o',vbr,pH,'-')
title('Titration of acetic acid')
xlabel('Vol NaOH [mL]')
ylabel('pH')
Finally we should plot the comparison. HOWEVER mine's not that wonderful. I'd hoped for better.
If you want a better fit, you need to "adjust" your constants. For example, Kw = 1.8e-14 gives a much better fit. (Note Kw is a function of temperature, and 1e-14 is not a golden constant! Similarly for your other constants.

类别

在 Help Center 和 File Exchange 中查找有关 Curve Fitting Toolbox 的更多信息

标签

产品


版本

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by