Fitting Data for the given equations below?

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

[John D'Errico]

So many questions.

Is pi 3.14159...? If not, it is a REALLY bad idea to use the variable pi, as pi is a rather useful number.

You talk about variables as being constants, but are they KNOWN constants? Or are they variables you must estimate?

You talk about theta as an unknown variable. It it just another unknown constant, to be estimated? What makes theta different from any other unknown constant, except that it appears in both equations? Is there an unknown value of theta for every data point?

You say you have data. How much data? How good is the data? What do you expect about the noise? Is there a good reason why you would not actually attach your data to the question? Is there any reason at all why not? (No, there cannot be any, at least if you want useful help. And if you don't want useful help, then why ask a question?) Do so as a .mat file, using the paper clip button to attach it to a comment or to your question. If you know the values of some of those constants, then tell us.

[Ashish Zore]

thanks for the information about pi. it is not 3.14 its a known variable called surface pressure (i.e. surface tension of solvent - surface tension of surfactant solution). this is a surface tension experiment. c is the concentration of the surfactant. And theta is the measure of amount adsorption of the surfactant at the air water interface (which causes surface tension reduction). Adsorption i.e theta cannot be measured (we dont have apparatus to measure it. it requires neutron scattering experiment and there are some technicalities with it). so the surface pressure pi and conc c are related by above two equations. the problem is dont know how to eliminate theta so that pi and c are related by substitution. so i am stuck with two equations. now i want to fit my data using these two equations using curve fitting and iteration.

Here is the data

c, mmol/l   pi, mN/m
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

Like i said pi (surface pressure) and c (concentration) are known variable and theta (surface adsorption) is unknown variable (it cannot be measured). the three are related as shown in the equations.

the unknown constants are a, b and omega. i want to determine those. R (gas constant) and T (temperature) are known constants.

to answer you question -What makes theta different from any other unknown constant, except that it appears in both equations? - theta is surface adsorption. it changes with concentration. higher the conc more is the adsorption. More the adsorption higher is the surface pressure pi. So theta is a variable i cannot measure.

About the .mat file, i have not yet analysed my data in matlab. i have pasted my data in this comment. I am sorry for not doing that.

the values of unknown constant omega and a can be guessed like a range omega - 0 to 1.0e+9 a = 0 to 20

the unknown constant b cannot be guessed.

Thank you. i appreciate you helping out. Please let me know if you want more information

[John D'Errico]

But what I wanted to get from you was what are the value of the unknown constants? What value does T have? How about R?

[Ashish Zore]

R = 8.314e+3 T=299K

Answer 1

[David Goodmanson]

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

Answer 2

[John D'Errico]

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.

You can download fminsearchbnd from the file exchange.

[John D'Errico]

Star, we have already been told:

 R (gas constant) and T (temperature) are known constants.

Which seems to make your comment a reiteration of known information. I can read. But what information are you adding? Unless you happen to know the values of those two constants, I'm sorry, but it appears to be nothing. Luckily, I did not need those values to solve most of the problem.

As for the actual value of R, yes, I could probably look it up. It is something that would depend on the units for this problem. But we don't know the temperature anyway.

[Ashish Zore]

 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.

Answer 3

[David Goodmanson]

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]

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.

Answer 4

[Ankit Kanthe]

Hi Ashish, Even I am trying to find the solution for the Frumkin adsorption isotherm. Did you get any success on that?

Thanks!

[Ashish Zore]

I was not able to figure out using matlab. Although I do have a software that solves the frumkin adsorption. Look up for the papers published by fairnerman, Miller, joos. They have also made a software to calculate for all the different isotherms.

[Ankit Kanthe]

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 : [e-mail] then that would be great.

Thanks for the help.

Answer 5

[Alex Sha]

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