Application Center - Maplesoft

App Preview:

Fitting vapor pressure data with regression

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application


 

Data Regression

> restart;

The problem at hand is the fitting of experimental data to equations used to represent the vapor pressure. We begin with a list of the temperatures (in Celcius) at which the vapor pressure has been measured.

> T(C) := [-36.7,-19.6,-11.5,-2.6,7.6,15.4,26.1,42.2,60.6,80.1];

[Maple Math]

Note the perhaps unconventional symbol used to denote the list of temperatures. As far as Maple is concerned [Maple Math] here is a name just like [Maple Math] would be, and can be assigned to any valid Maple object; in this case a list of temperature values. The corresponding values of the vapor pressure (in mm Hg) is next

> P(mmHg) := [1,5,10,20,40,60,100,200,400,760];

[Maple Math]

A quick check to see that the number of values of temperature and pressure are equal.

> nops(P(mmHg)), nops(T(C));

[Maple Math]

We plot the data to see its general form

> datapoints := [seq([T(C)[k],P(mmHg)[k]],k=1..nops(T(C)))];

[Maple Math]

> dataplot:=plot(datapoints,T=-60..100,style=point,symbol=box): dataplot;

[Maple Plot]

We continue by attempting to fit simple polynomials to the data. We define a function to give us a polynomial in [Maple Math] of any order.

> p := n->sum(a[k]*T^k,k=0..n);

[Maple Math]

A quadratic, for example, is given by

> p(2);

[Maple Math]

and a cubic is

> p(3);

[Maple Math]

Maple includes a least squares fitting procedure in its statistical package. The package is loaded as follows

> with(stats):

We fit a quadratic to the data as follows

> n:=2;
expr := P=p(n);

[Maple Math]

[Maple Math]

> eqn[n]:=fit[leastsquare[[T,P], expr, {seq(a[k],k=0..n)}]]([T(C),P(mmHg)]);

[Maple Math]

We can fit higher order polynomials just as easily

> n:=3; expr := P=p(n);

[Maple Math]

[Maple Math]

> eqn[n]:=fit[leastsquare[[T,P], expr, {seq(a[k],k=0..n)}]]([T(C),P(mmHg)]);

[Maple Math]

> n:=4; expr := P=p(n);

[Maple Math]

[Maple Math]

> eqn[n]:=fit[leastsquare[[T,P], expr, {seq(a[k],k=0..n)}]]([T(C),P(mmHg)]);

[Maple Math]

> n:=5; expr := P=p(n);

[Maple Math]

[Maple Math]

> eqn[n]:=fit[leastsquare[[T,P], expr, {seq(a[k],k=0..n)}]]([T(C),P(mmHg)]);

[Maple Math]

We can compare these expressions graphically

> n:='n':
polyplot:=plot([seq(rhs(eqn[n]),n=2..5)],T=-40..100,color=[red,black,blue,green]):

> plots[display]({polyplot,dataplot});

[Maple Plot]

From which we observe that the red line (the quadratic) is clearly poorer than the others; all of which seem more or less equally good. However, if we plot the equations on a semilog plot we see that the polynomials suffer from a rather serious problem.

> lplot1:=plots[logplot]([seq(rhs(eqn[n]),n=2..5)],T=-40..100,color=[red,black,blue,green]):

> lplot2:=plots[logplot](datapoints,style=point,symbol=box):

> plots[display]({lplot1,lplot2});

[Maple Plot]

From this diagram (if viewed in color) we would notice that the fourth and fifth order polyomials are best (as might be expected).

We now attempt to fit the Clausius Clapeyron (CC) equation to the data.

> CCeqn:= log10(P) = A-B/(T): CCeqn;

[Maple Math]

The bottom line of the second term on the right is the temperature in Kelvin. We define two new variables

> v1:= Tau =1/T: v1;

[Maple Math]

> v2:=rho=log10(P): v2;

[Maple Math]

and substitute into the CC equation.

> CC2:=subs(rhs(v2)=lhs(v2),rhs(v1)=lhs(v1),CCeqn): CC2;

[Maple Math]

Thus, this equation is linear in [Maple Math] and [Maple Math] .

The next step is to transform the data, beginning with the vapor pressures

> lnpdata := evalf(map(log10,P(mmHg)));

[Maple Math]

and now for the temperatures.

> ToKelvin:=x->x+273.15;
T(K) := map(ToKelvin,T(C));

[Maple Math]

[Maple Math]

We must take the reciprocal of these values

> r := x-> 1/x;
Tover := map(r,T(K));

[Maple Math]

[Maple Math]

We are now ready to fit the transformed data

> CCfit:=fit[leastsquare[[Tau,rho], CC2, {A,B}]]([Tover,lnpdata]): CCfit;

[Maple Math]

This result can be expressed in terms of our original variables

> CC3:=subs(v1,v2,T=ToKelvin(T),CCfit): CC3;

[Maple Math]

We make a function of the right hand side for plotting purposes

> CCF:=unapply(10^rhs(CC3),T);

[Maple Math]

and display this new result along with the data.

> n:='n':
p1:=plots[logplot]([CCF(T)],T=-40..100,color=[red]):

> plots[display]({p1,lplot2});

[Maple Plot]

From this we see immediately that the Clausius Clapeyron form is far superior to the simple polynomial fits.

Another popular equation used for fitting vapor pressure data is the Antoine equation.

> Antoine := log10(P)=A-B/(T+C): Antoine;

[Maple Math]

This equation cannot be made linear in [Maple Math] and [Maple Math] . However, some rearrangement of the above expression gives:

> A1:=collect(expand(normal(Antoine*(T+C))/T),T):

> t:=op([1,2],A1):

> A2:=A1-(t=t): A2;

[Maple Math]

This equation is linear in [Maple Math] and [Maple Math] with parameters [Maple Math] , [Maple Math] , and [Maple Math] .

We define new temperature and pressure variables as was done for the Clausius-Clapeyron equation.

> v1:=Tau=1/T: v1;

[Maple Math]

> v2:=rho=log10(P): v2;

[Maple Math]

We must also define the combined variable:

> v3:=Gamma = rho*Tau: v3;

[Maple Math]

With the above definitions we may rewrite the Antoine equation as

> read `e:/maple/utils/utils.mpl`:

> A3:=subs(Solve([v1,v2,v3],[T,ln(P),Tau]),A2):
A4:=algsubs(Tau*rho=Gamma,A3):
A4:=sort(A4,[a[0],a[1],a[2]]): A4;

[Maple Math]

We define new parameters:

> params := [a[0]=A,a[1]=coeff(rhs(A4),Tau,1),a[2]=-coeff(rhs(A4),Gamma,1)]: op(params);

[Maple Math]

and rewrite the equation in terms of these parameters:

> A4:=subs(map(flip,params,equation),A4): A4;

[Maple Math]

The next step is to obtain the necessary data. We already have the column of [Maple Math] ; the column of values of [Maple Math] is obtained as follows.

> r := x-> 1/x;
Tover := map(r,T(K));

[Maple Math]

[Maple Math]

> Gammavalues := [seq(Tover[k]*lnpdata[k],k=1..nops(lnpdata))];

[Maple Math]

We are now ready to fit the transformed data

> Antoinefit:=evalf(fit[leastsquare[[Tau,Gamma,rho], A4, {a[0],a[1],a[2]}]]([Tover,Gammavalues,lnpdata])):
Antoinefit;

[Maple Math]

In the next few lines we reconstruct the Antoine equation in terms of oroginal variables:

> peqn1:=a[0]=coeff(coeff(rhs(Antoinefit),Tau,0),Gamma,0): peqn1;

[Maple Math]

> peqn2:=a[1]=coeff(rhs(Antoinefit),Tau): peqn2;

[Maple Math]

> peqn3:=a[2]=-coeff(rhs(Antoinefit),Gamma): peqn3;

[Maple Math]

> params;

[Maple Math]

> params2:=solve(subs(peqn1,peqn2,peqn3,{op(params)}),{A,B,C}): params2;

[Maple Math]

This result can be expressed in terms of our original variables

> A5:=subs(params2,T=ToKelvin(T),Antoine): A5;

[Maple Math]

We make a function of the right hand side for plotting purposes

> AF:=unapply(10^rhs(A5),T);

[Maple Math]

and display this new result along with the Clausius-Clapeyron equation and the data.

> p1:=plots[logplot]([CCF(T),AF(T)],T=-40..100,color=[red,blue]):

> plots[display]({p1,lplot2});

[Maple Plot]

In this case the Antoine fit (in blue) appears to be even better than the was obtained with the Clausius-Clapeyron equation (in red). Contrast the behavior of both methods with that of the simple polynomial fits given above.

>