MATLAB Answers

0

Fitting model to titration data

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 Comments

Sign in to comment.

Products


Release

R2019a

2 Answers

Answer by Rafael Ibáñez on 13 Apr 2019
 Accepted Answer

Thank you very much

  0 Comments

Sign in to comment.


Answer by David Wilson on 13 Apr 2019
Edited by David Wilson on 13 Apr 2019

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.

  1 Comment

Thnaks¡¡¡¡
It works

Sign in to comment.