# Fitting Data for the given equations below?

8 views (last 30 days)
Ashish Zore on 28 Jul 2017
Commented: Ankit Kanthe on 28 Feb 2018 I am fitting my experimental data with the equations above. In the equation, pi, theta and c are variables and a, b and omega are constants. However, of the three variable only pi and c are known as i have measured them experimentally. theta is an unknown variable. I want to fit my data to these equations and get the values of the constants. i cannot eliminate the unknown variable theta by substitution. Rather i dont know how to do that. i want to know how to fit my data using matlab with those equations. its a bit complicated than usual equations.
thank you
Ashish Zore on 31 Jul 2017
R = 8.314e+3 T=299K

David Goodmanson on 29 Jul 2017
Edited: David Goodmanson on 31 Jul 2017
Hello Ashish,
SEE SECOND ANSWER ALSO POSTED I wanted to revise this one but was having technical difficulties
I am going to take at face value the idea that theta is unknown, and that one way or another all variables that are not theta are known. Yes it would be a bad idea to actually use 'pi' as a variable, it should be renamed, but for the moment, defining:
A = (pi*w)/(R*T)
B = bc
then
-A = log(1-th) + a*th^2
B = (th/(1-th))*exp(-2*a*th )
and this is to be solved for theta. Exponentiating the first equation and multiplying by the second leads to
B*exp(-A) = th.*exp(-2*a*th +a*th.^2)
with the solution
a = 6; % for example
th = .7;
% make some variables with a known outcome
A = -( log(1-th) + a*th^2 )
B = (th/(1-th))*exp(-2*a*th )
% solve
fun = @(th) th.*exp(-2*a*th +a*th.^2) - B*exp(-A);
th = fzero(@(th) fun(th), .5)
% check for smallness
A + log(1-th) + a*th^2
B - (th/(1-th))*exp(-2*a*th)
th = 0.7000
diff1 = -1.7764e-15
diff2 = -6.5052e-19

John D'Errico on 29 Jul 2017
Edited: John D'Errico on 29 Jul 2017
Again, pi is a BAD name for a variable. I won't use it. I'll call it P instead.
cp = [4.09 19.66
3.57 16.66
3.05 15.06
2.54 13.46
2.02 11.06
1.77 9.96
1.51 8.66
1.26 7.46
1.01 6.46
0.75 4.26
0.50 3.26
0.25 1.36];
c = cp(:,1);
P = cp(:,2);
I have not been given the value of R or T. So I can't solve the problem completely. But omega is an unknown constant too. So I can freely define a new variable as
K = R*T/omega
I'll estimate that ratio instead. That leaves us two equations for each data point:
P = K*(log(1-theta) + a*theta^2)
b*c = theta/(1-theta)*exp(-2*a*theta)
Note that the natural log function in MATLAB is log, NOT ln. Take the log of the second equation.
log(b) + log(c) = log(theta) - log(1-theta) - 2*a*theta
Why do I take the log here? Because it is terribly important to do so! If a is at all large, then exp(-2*a*theta) will cause an exponent underflow. We will get a zero out, and therefore garbage for a result.
Next, we can write the first equation as
P/K = log(1-theta) + a*theta^2
Ok, so we can take it that far pretty easily. Now, lets first check one thing. Do we have sufficient information to infer the values of all variables?
We have 12 data points, thus 12 combinations of c and P. For each combination of c and P, theta is unknown. So there are 12 unknown values of theta. As well, we don't know the values of a and b, nor do we have the value of K. (Of course, once you give me the values of R and T, I could infer the value of omega.)
So for 12 data points, we have 12+3=15 unknown parameters, and 12*2=24 equations that relate these parameters together.
So in theory, there is sufficient information to estimate all 15 unknowns. Of course, theory and practice often seem to get in the way of each other.
What do we know about the unknown parameters?
We know that theta lies in the interval (0,1). It cannot lie outside of that interval, since otherwise we will see the log of zero or a negative number.
a you said should lie in the interval [0,20].
You told us a range for omega, but since you gave us no information about R or T, we are lost there.
b is also a complete unknown. You gave nothing at all there, except we know it must be positive.
I think my approach will be to use an optimization tool to pose various values for the unknown constants a, b, and K (and therefore implicitly omega.) Given any value of those three parameters, I now have TWO equations for each value of theta. Find the value of theta in the interval (0,1) which minimizes the errors for those two equations, for EVERY pair (c,P). Once that is done, return the overall sum of squares of residuals to the optimization tool.
This scheme effectively optimizes all 15 parameters at once.
Now, let me take a stab at it...
function [a,b,K,theta] = optimizeparameters(c,P)
% Lacking intelligent starting values for those parameters,
% I'll just pick some.
abk0 = [10,1,1];
% Now call the optimizer. I'll just use fminsearchbnd,
% this will allow me to constrain a, b, and K to be positive.
% I won't set any upper bounds.
LB = [0 0 0];
UB = [inf inf inf];
options= optimset('fminsearch');
options.Display = 'iter';
abkfinal = fminsearchbnd(@abkfun,abk0,LB,UB,options);
% once we know the final values of a, b and K, get theta for those values.
[~,a,b,K,theta] = abkfun(abkfinal);
% a nested function to compute a sum of squares of residuals,
% as well as theta.
function [ssq,a,b,K,theta] = abkfun(abk)
a = abk(1);
b = abk(2);
K = abk(3);
n = numel(c);
% solve for each theta.
theta = zeros(size(c));
ssq = 0;
for i = 1:n
thetafun = @(theta) (log(1-theta) + a*theta^2 - P(i)/K)^2 + ...
(log(b) + log(c(i)) - log(theta) + log(1-theta) + 2*a*theta).^2;
[theta(i),ssqi] = fminbnd(thetafun,eps,1-eps);
ssq = ssq + ssqi;
end
end
end
Now, lets run the code passing in c and P.
[a,b,K,theta] = optimizeparameters(c,P)
(I've cut away the fminsearch display output here...)
a =
3095.1
b =
8.8421e-11
K =
664.46
theta =
0.0025503
0.0025737
0.0026008
0.0026322
0.0026713
0.0026936
0.0027204
0.0027509
0.0027882
0.0028384
0.0029067
0.0030234
I'll note that you said a should be up to about 20. But I question if you know that limit, since a seems to want to go to around 3000.
Given the value of K, IF you know R and T, then omega is trivial to recover.
##### 2 CommentsShowHide 1 older comment
Ashish Zore on 31 Jul 2017
R = 8.314e+3 T=299K
range for omega = 0 to 1.0e+9,
i made a mistake about a, it can also be negetive so ,a = -20 to 20. this is from other experimental data.

David Goodmanson on 31 Jul 2017
Edited: David Goodmanson on 31 Jul 2017
Hello Ashish
[ I wanted to append to my previous answer but for some reason the site would not let me copy blocks of text into that answer ]
Now that I know the ground rules (I think), here is another attempt, although it only fits two of three parameters. I assume that [a,b,w] are unknown constants to be fitted; theta is also an unknown quantity; that you know R and T, and you have a curve of p(Pi) vs. c. Correct?
T is known so I arbitrarily assumed T = 294 here. There are two equations, four unknowns. If you could eliminate theta you would have one equation, three unknowns, so the only way to get those is to pull them out of the behavior of p vs. c. But if you look at the curve, it’s pretty smooth (good data from that point of view), and it does not look like you can accurately fit more than two parameters. The code
% mmol/l mN/m
x = ...
[4.09 19.66
3.57 16.66
3.05 15.06
2.54 13.46
2.02 11.06
1.77 9.96
1.51 8.66
1.26 7.46
1.01 6.46
0.75 4.26
0.50 3.26
0.25 1.36]
c = x(:,1); % mol/m^3
p = x(:,2)*1e-3; % N/m
a = polyfit(p,c,2)
fig(1)
plot(p,c,p,polyval(a,p))
shows a very nice fit to a parabola that passes through the origin. (I changed to m^3 since I’m a physics guy, not a chemist). With just the two parameters I decided (maybe you have done this already) to try a simpler problem: assume a = 0. Then you can eliminate theta, leading to
c/c0 = exp(p w /RT) -1
where the concentration c0 is the inverse of your constant b. This could be fitted in a fancy way but the curve of p vs c is a pretty shallow parabola so the exponent can be expanded out, leading to
c/c0 = (p w /RT) + (p w /RT)^2/2.
This says the curve should pass through the origin and be concave upward which it is. With the polynomial fit and a bit of algebra the solution is shown below.
RT = 8.31*294;
a = polyfit(p,c,2)
d2 = a(1)
d1 = a(2)
w = 2*(d2/d1)*RT
c0 = (1/2)*(d1^2/d2)
%results
w = 1.2356e+05 % m^2/mol
c0 = 2.8193 % mmol/l
For w, the intermolecular spacing on the surface would be
s = sqrt(w/6e23) s = 4.5379e-10
or 4.5 Angstroms which seems reasonable. I had no idea what c0 might be but it is comparable to your measured values for c. The plot of theta = c/(c + c0) varies from .1 to .6.
Ashish Zore on 31 Jul 2017
I have already done with assumption of a=0. that is a simpler model called langmuir model. it give a single equation since theta can be eliminated. what i am trying to solve is Frumkin model where a is not zero.

Ankit Kanthe on 27 Feb 2018
Edited: Ankit Kanthe on 27 Feb 2018
Hi Ashish, Even I am trying to find the solution for the Frumkin adsorption isotherm. Did you get any success on that?
Thanks!
##### 2 CommentsShowHide 1 older comment
Ankit Kanthe on 28 Feb 2018
Thanks for your reply. Yes I am following the Miller papers. If possible can you forward me the software link or the article to my email id : kantheankit@gmail.com then that would be great.
Thanks for the help.

Alex Sha on 11 Oct 2019
seems too complicated for Matlab, the follows are the code of 1stOpt, another software package:
Constant RT = 8.314E+3*299;
ParVariable theat;
SharedModel;
Variable p,c;
Function p = -RT/omega*(ln(1-theat)+a*theat^2);
c = theat/(1-theat)*exp(-2*a*theat)/b;
Data;
4.09 19.66
3.57 16.66
3.05 15.06
2.54 13.46
2.02 11.06
1.77 9.96
1.51 8.66
1.26 7.46
1.01 6.46
0.75 4.26
0.50 3.26
0.25 1.36
Results:
Parameter Best Estimate
-------------------- -------------
omega -492023.278565147
a -0.0356831451208793
b -0.0291587333259045
theat0 -5.42701040750307
theat1 -1.10987296581655
theat2 -0.878195386477076
theat3 -0.701158079125353
theat4 -0.502244639937896
theat5 -0.427446952224332
theat6 -0.349527045443767
theat7 -0.285405918079108
theat8 -0.236526830070603
theat9 -0.143995617121495
theat10 -0.105885491091172
theat11 -0.0416372367071209  