Fitting model to titration data
9 次查看(过去 30 天)
显示 更早的评论
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
0 个评论
采纳的回答
更多回答(1 个)
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. ![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/213597/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/213597/image.png)
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.
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!