Numerical Method for Analysis of Surface ... - ACS Publications

3630

Langmuir 1996, 12, 3630-3642

Numerical Method for Analysis of Surface Heterogeneity in a Case of Finite Diversity of Adsorption Sites Vadim Sh. Mamleev* and Esen A. Bekturov A. B. Bekturov Institute of Chemical Sciences, Kazakh Academy of Sciences, Valikhanov str. 106, Almaty 480100, Kazakhstan Received March 6, 1995. In Final Form: January 2, 1996X A new algorithm is suggested for the calculation of energies and fractions of adsorption sites of different kinds from adsorption isotherms. The method is recommended for the investigation of crystalline surfaces when the number of distinguishable adsorption sites is not too large (=1-5). The algorithm was applied in the treatment of the adsorption isotherm of nitrogen on the surface of maximally hydroxylated silica at 78 K. The procedure, which allows one to calculate the isosteric heat of adsorption with a correction for effects of multilayer adsorption, is used for calculation of the energies of adsorption.

1. Introduction The methods of quantitative estimation of the energetic heterogeneity of surfaces were often severely criticized. The history of the problem and main physical postulates accepted in the theory of heterogeneity have been considered in a number of surveys.1-3 It is necessary to solve, as a rule, ill-posed (or “incorrect”4,5) problems to study heterogeneity of adsorbents from adsorption isotherms. Earlier6 we considered an algorithm for the heterogeneity analysis in a case of continuous distributions of adsorption sites (centers) with respect to their adsorption energies. The purpose of the present work is to discuss an algorithm that can be applied to the case of a discrete distribution. We shall imply below that any energetic spectrum consisting of a sum of peaks with widths =RT (or less) can be considered as a discrete distribution because the representation of a narrow peak in form of the deltafunction does not lead to an appreciable accuracy loss. For the quantitative description of the adsorption on heterogeneous surfaces, two limiting assumptions are usually taken. The first one is that the surface consists of a multitude of uniform “patches” (the patchwise model). If interaction between molecules adsorbed on different “patches” may be neglected, then the experimental adsorption isotherm is represented as a sum of local isotherms corresponding to individual “patches”. The second assumption is that the adsorption sites on the surface are randomly “mixed” (the model of random heterogeneity). In the latter case an overall isotherm is also represented as a sum of the local ones but includes a correction for nonadditivity. The models of patchwise and random heterogeneity may be written7 as N

Fiθ(P,i) ) n(P) ∑ i)1

na

(1)

N

Fiθ(P,i, Θ) ) n(P) ∑ i)1

na

(2)

where n(P) is the adsorbate quantity per unit of surface * To whom correspondence should be addressed. FAX: 7-327261-57-65. E-mail: [email protected] X Abstract published in Advance ACS Abstracts, May 15, 1996. (1) House, W. A. Colloid Sci. 1983, 4, 1. (2) Jaroniec, M.; Braüer, P. Surf. Sci. Rep. 1986, 6, 65. (3) Sircar, S.; Myers, A. L. Surf. Sci. 1988, 205, 353. (4) Tikhonov, A. N.; Arsenin, V. Ya. The Methods of Solution of Incorrect Problems; Nauka: Moscow, 1979 (in Russian).

S0743-7463(95)00181-8 CCC: $12.00

area, na is an adsorption capacity, N is the number of the kinds of the distinguishable adsorption sites on the surface, i and Fi are the adsorption energy and the fraction of the adsorption sites belonging to the kind with the number i, θ is the degree of occupation of these sites (the local isotherm), P is a partial pressure of adsorbate in a bulk phase, Θ(P) ) n(P)/na is the overall degree of occupation of the monolayer by adsorbate molecules. The difference between the models with patchwise and random heterogeneity, in the quantitative form probably for the first time, has been investigated by Hill.7 Equations 1 and 2 derived by Hill for localized adsorption were commonly used in different publications1-3 including present-day ones.8-10 It has not been until recently11 that assumptions of Hill’s theory have been called in question. The principle of additivity (see eq 1) for the patchwise surfaces, e.g., for the surfaces of monocrystals, seems obvious. However, eq 2 raises doubts because it neglects the correlation between pair distribution of nearest adsorption sites of different kinds and spatial distribution of interacting pairs of adatoms. In the mean field theory (Bragg-Williams aproximation) this correlation is ignored, and as natural consequence, eq 2 is obtained. When using, e.g., quasichemical approximation in treating the local isotherm, one can obtain other equations. However, the deductions of Hill7 probably contain an inaccuracy. Assuming that adsorption sites are randomly “mixed”, he has obtained eq 2 not only by employing the mean field theory but also by employing the quasichemical approximation. By using more strict statistical methods in the quasichemical approximation for a lattice model, Tovbin has derived11 more correct equations. Note that7 eqs 1 and 2 and slightly more accurate equations derived until now11 are far from being perfect. At the same time the use of more complicated theory would not lead to discovery of any new qualitative effects. Therefore the use of simple equations seems reasonable. We deem that the mean field theory, at present, optimally solves the compromise between simplicity and accuracy. (5) Tikhonov, A. N.; Goncharskii, A. V.; Stepanov, V. V.; Yagola, A. G. Regularization Algorithms and a priori Information; Nauka: Moscow, 1983 (in Russian). (6) Mamleev, V. Sh.; Bekturov, E. A. Langmuir 1996, 12, 441. (7) Hill, T. L. J. Chem. Phys. 1949, 17, 762. (8) Koopal, L. K.; Vos, C. H. W. Langmuir 1993, 9, 2593. (9) Szombathely, M. v.; Braüer, P.; Jaroniec, M. J. Comput. Chem. 1992, 13, 17. (10) LamWan, J. A.; White, L. R. J. Chem. Soc., Faraday Trans. 1991, 87, 3051. (11) Tovbin, Yu. K. Theory of Physical Chemistry Processes at GasSolid Interface; Nauka: Moscow, 1990 (in Russian).

© 1996 American Chemical Society

Analysis of Surface Heterogeneity

Langmuir, Vol. 12, No. 15, 1996 3631

We shall suppose that the function n(P) (an overall isotherm) is measured in the physical experiment. As to the local isotherm θ(P,) (or θ(P,,Θ)), it is chosen on the basis of proper physical models of adsorption. Such information is sufficient to calculate from eq 1 (or eq 2) the two sets of the values {Fi} and {i}. If nR is known, one may divide both parts of eq 1 (or eq 2) by na to find Θ(P). Otherwise one should search for the values {naFi}. Taking account of the normalization N

Fi ) 1 ∑ i)1 na is calculated by summing up N

na )

∑ i)1

N

Fi ∑ i)1

naFi ) na

To provide for the normalization it is sufficient to divide naFi (1 e i e N) by na. In the case of continuous distribution, sums 1 and 2 are replaced by the integrals

∫ F() θ(P,) d ) n(P)

(3)

∫ F() θ(P,,Θ) d ) n(P)

(4)

b

na

a

b

na

a

where F() is the statistical density of distribution of adsorption sites in accordance with adsorption energy; , a, and b are the adsorption energy and its minimal and maximal values. Much work1-3 has been devoted to the problem of the solution of the integral equations 3 and 4. In particular, algorithms6,8,9 have been developed where regularization is used for reproducibility of calculated F() and the necessary constraint

(a e e b)

F() g 0

(5)

of nonnegativity of this function is taken into account for physical correctness of solution. It is obvious that in a case of discrete distribution instead of constraint 5, we have

Fi g 0

(1 e i e N)

If F() is the sum of Dirac-delta functions N

F() )

Fiδ( - i) ∑ i)1

eqs 3 and 4 turn into eqs 1 and 2. Consequently, eqs 3 and 4 have a more general form in comparison to eqs 1 and 2. Nevertheless the algorithm of the calculation of discrete distribution is sure to be of value because one may foresee the situations when representations 1 and 2 are preferable in comparison with eqs 3 and 4. This question needs to be discussed in brief. Let us represent eq 3 by using the integral of Stieltijes

∫ F() θ(P,) d ) na∫01θ(P,(F)) dF

n(P) ) na

b

a

(6)

where F is the integral adsorption energy distribution

F)

∫F() d a

Figure 1. Results of evaluating the distribution consisting of the two narrow peaks (a, b) which one may expect when using small (a, a′) and large (b, b′) numbers of the equidistant nodes and when using the algorithm of calculation of discrete distributions (c, c′) with N ) 2. The use of the integral distribution functions (a′-c′) is a common method to display both discrete and continuous spectra.

Any integral 6 can be approximated as a finite sum using, e.g., the rule of rectangles N

∫01θ(P,(F)) dF ≈ na∑θ(P,i) ∆Fi

n(P) ) na

(7)

i)1

where {i} are the nodes on the energy axis, ∆Fi ) Fi Fi-1, (1 e i e N), F0 ) F(R) ) 0. Having calculated the set {∆Fi}, one may find the knot values of F() by the formula

F(i*) ≈ (Fi - Fi-1)/(i - i-1)

(8)

where i* ) (i + i-1)/2. Note that eq 7 turns into eq 1 if denoting Fi ) ∆Fi. Thus, the problems of calculation {i} and {Fi} or {i} and {∆Fi} are equivalent. The first question to be solved is how many nodes should be chosen on the energy axis to provide for the necessary accuracy of the numerical integration by eq 7. Let us assume that distribution to be found represents two narrow peaks having the widths =RT (or Tc will be considered below. In this case both local isotherms 10 and 11 are continuous, monotonic functions of pressure. It is important to generalize the method of heterogeneity analysis to include multilayer adsorption. Unfortunately, no satisfactory theory of multilayer adsorption is available even for the uniform surfaces. The theory of Brunauer, Emmett, and Teller (BET) is used in practice most often. Nobody has proved correctness of the BET theory for uniform surfaces. Moreover, good reasons exist to assert categorically that this theory is not suited for description of multilayer adsorption even in the simplest systems such as noble gases on a uniform surface of carbon black.20 Thus, the BET equation cannot be used as a local isotherm. Though in a number of publications it has been used as the kernel of the integral equation (see surveys1,2). At the same time it is known that the BET theory approximates well the overall isotherms of adsorption on heterogeneous surfaces; i.e., the BET equation can be used as the integral isotherm. Following this theory one can introduce the simplest correction21

ΘB(P) ) Θ(P)/(1 - P/P0)

(14)

or

Θ(P) ) ΘB(P)(1 - P/P0) where ΘB(P) is an experimental isotherm of multilayer adsorption, Θ(P) is an integral isotherm of monolayer adsorption, and P0 is a pressure of saturated vapors of (19) Jaycock, M. J.; Parfitt, G. D. Chemistry of Interfaces; J. Wiley & Sons: New York, 1981. (20) Thomy, A.; Duval, X.; Regnier, J. Surf. Sci. Rep. 1981, 1, 1. (21) Osćik, J. Adsorption; Polish Scientific Publishers. PWN: Warsaw, 1982.

adsorbate at the temperature T < T0. Note that formula 14 is formally derived if the local BET isotherm is employed as the kernel but the constant C for all sites is large enough (C . 1). Equation 14, of course, must be considered as an empirical dependence. Nevertheless at this stage of the theory development it is difficult to choose a more suitable equation for recalculation of ΘB(P) into Θ(P). At least it gives a correct limit of the isosteric heat qst:qst f L as P f P0, where L is the heat of vaporization of liquid adsorbate. It is worthwhile to note that the present work is not intended for development of any new thermodynamic concepts. As the theory will advance, other isotherms can be used for the analysis. 2. The Numerical Method Let us assume that the values of measured function n(P) are known in M points: n(Pm)(1 e m e M). If n(P) is measured within the experimental interval of pressure, P1 e P e PM with, approximately, the constant relative error and the number of measurements M is large enough, then the optimal sets {Fi}, {i} must ensure the minimum of the following function of 2 × N variables: M

∆[{Fi},{i}] ) M

-1

N

{1 - [na/n(Pm)]∑Fiθ(Pm,i)}2 ∑ m)1 i)1

The value of ∆1/2 is the mean relative deviation (MRD) between theoretical and experimental integral isotherms. It is obvious that we have to solve the problem

∆[{Fi},{i}] w min Fi g 0

(1 e i e N)

(15) (16)

Let us introduce the new function

Gi(P) ) Fiθ(P,i) or Gi(P) ) Fiθ(P,i,Θ) to rewrite eqs 1 and 2 in the form N

Gi(P) ) n(P)/na ∑ i)1

(17)

The members of series 17, i.e., the functions Gi, can be numbered by arbitrary manner. Thereby the number of the minima of ∆[{Fi},{i}] is divisible by N!. The physical uniqueness of a solution denotes that the number of the minima is equal to N! and all minima have one and the same “depth”. Below, we shall imply that such minima are actually one minimum. With an integral isotherm being set exactly, by using the method of the Pade´ approximation, Oh and Kim17 have proved the uniqueness of solution in a case of the Langmuir local isotherm. However, uniqueness may be violated if an integral isotherm includes errors. With the local isotherms 10 and 11 the solution, probably, can lose uniqueness even if an integral isotherm is known exactly. Although the solution is likely to be unique in all cases if N is not too large (N = 1-5). The necessary (but insufficient) condition of the minimum 15, 16 consists in existence of such nonnegative (λi

3634

Langmuir, Vol. 12, No. 15, 1996

Mamleev and Bekturov

g 0) Lagrange multipliers which provide for the following system:

∂∆/∂Fi - λi ) 0 λiFi ) 0 ∂∆/∂i ) 0

(1 e i e N) (1 e i e N) (1 e i e N)

(18) (19) (20)

Equations 18-20 are often named22 the Kuhn-Tacker conditions. Let the sets {Fi}, {i} be numbered by such a way that, as a result of minimization, the first J components {Fi} (1 e i e J) would be positive (the inactive22 constraints 19) and the remainder of the components {Fi} (J < i e N) would be equal to zero (the active22 constraints 19). It is seen that to find a physically correct solution, it is enough to abbreviate the number of the sought for components to J. Then, according to conditions 19, the multipliers λi will be equal to zero when 1 e i e J and the function ∆ will have in the region Fi > 0 (1 e i e J) the absolute (unconditional) minimum that can be found by solving an algebraic system of equations satisfying the equalities

∂∆/∂Fi ) 0,

∂∆/∂i ) 0 (1 e i e J)

realization of the latter method the positive term ensuring the fulfillment of condition 16 is added to function ∆

(21)

If the model were exact and experimental values of n(P) were not to include the errors, then by increasing the number N in sequence 1, 2, 3, ..., one might find the critical number N ) J for which a solution of the system 21 corresponds to zero minimum of ∆. If the minimum of ∆ is equal to zero and it is reached with the positive components {Fi}, there is no necessity to envisage constraint 16 in the algorithm. However, the input information is inaccurate. It is well-known that because of this inaccuracy the amplitudes {Fi} at the zero-minimum point may oscillate between large positive and large negative values.9,23 Thus, in the absence of a priori physical information about the number N, constraint 16 surely must be taken into account. In a number of publications8,9,15 for taking into account the nonnegativity of solution of the integral equations 3 and 4, the NNLS (nonnegative least squares) method of Lawson and Hanson24 has been used. Some circumstances hamper the application of this method for solving the problem in question. First, the NNLS method is destined for finding only a quadratic minimum. Approximate, reiterative replacement of ∆ by a quadratic function and a subsequent minimization of this approximate function can take away from the minimum of the exact function rather than get nearer to this minimum. Second, the NNLS algorithm seems to be not able to transfer the active constraints into the set of inactive ones. Note that if some component Fi has been drawn into the set of active constraints, then to transfer it into a set of inactive constraints one must calculate the partial derivative ∂∆/ ∂i. However, this derivative will always be zero if Fi ) 0. In other words, the positions of the points with zerovalued amplitudes do not influence the value of ∆ and can be arbitrary, therefore, even with known multipliers {λi} eqs 18-20 do not allow one to find the energies i (J < i e N). We find it hard to propose any modification of the NNLS algorithm to overcome the above difficulties. At the same time the method of barrier functions6,22 seems to be suited ideally for solving the present problem. For (22) Gill, P. E.; Murray, W.; Wright, M. H. Practical Optimization; Academic Press: London, 1981. (23) Merz, P. J. Comput. Phys. 1980, 38, 64. (24) Lawson, C. L.; Hanson, R. J. Solving Least Squares Problems; Prentice-Hall: Englewood Cliffs, NJ, 1974.

N

∆Rs[{Fi},{i}] ) ∆[{Fi},{i}] + Rs

[ln(Fi) + D]2 ∑ i)1

(22)

where Rs is a small positive parameter and D is a positive number (D = 10). When Rs f 0 the minima of ∆Rs[{Fi},{i}] and ∆[{Fi},{i}] coincide (see the Appendix); i.e., when Rs f 0 the conditional minimization of the function ∆[{Fi},{i}] subject to constraint 16 is actually equivalent to the absolute minimization of ∆Rs[{Fi},{i}], viz., the problem 15, 16 comes to

∆Rs[{Fi},{i}] w min

(23)

or

∂∆Rs/∂Fi ) 0, ∂∆Rs/∂i ) 0 (1 e i e N) Let {Fhi},{ji} and {F˜i}, {˜i} be the sets of {Fi},{i} which correspond to the minima 15, 16 and 23 respectively. One may expect that F˜i ≈ Fhi for inactive constraints and F˜i = 0 for active constraints. When choosing Rs, it is natural to proceed from the requirement that perturbations caused by the barrier additive would be small. For example, Rs can be taken so that the quantities |∂∆/∂Fi| (1 e i e J) corresponding to the inactive constraints would not exceed the chosen small parameter δR. Taking account of equalities

∂∆Rs/∂Fi ) ∂∆/∂Fi + 2Rs[ln(F˜i) + D]/F˜i ) 0 we have the estimation

Rs = 0.5|δRFhj/[ln(Fhj) + D]| where j is the index of the component with the minimal derivative |∂∆/∂Fj| ) δR. When Rs f 0 we have δR f 0. In our calculations (see the fourth section) we used either Rs ) 10-6 (δR = 10-4-10-3) or Rs ) 0. By comparing the derivatives of ∆ and ∆Rs, one may estimate the ordinates corresponding to the active constraints

2Rs[ln(F˜i) + D]/F˜i ≈ -λi or

F˜i ≈ exp(-F˜iλi/2Rs) exp(-D) Thus, when Rs f 0, instead of zero-valued ordinates we have the negligibly small ordinates F˜i , exp(-D) (J < i e N). In spite of this smallness, nevertheless, the corresponding derivatives ∂∆/∂i (J < i e N) are not equal to zero, as a consequence, some energies ˜ i (J < i e N) can result from minimization 23, though, if F˜i , exp(-D) it is not worthwhile, of course, to ascribe to the energy ˜ i any physical sense. It is noteworthy that occurrence of active constraints actually signifies a degeneracy of distribution. Indeed, existence of any two components with Fhi*> 0, ji and Fhj ) 0, arbitrary jj is equivalent to existence of one component with the ordinate Fhi + Fh*j and the energy ji. Such a situation is displayed in Figure 3.


Langmuir, Vol. 12, No. 15, 1996 3635

The criterion of degeneracy can be written as

∆[{Fhi},{ji}] - ∆[{Fhi},{ji},FhN+1,jN+1] , ∆ h

Figure 3. Result of approximation (the solid line (a)) of the integral isotherm calculated by formula 24 (points (a)) in a situation when the artificial distribution (b) has one positive and one negative ordinate: c, N ) 1; d, e, N ) 2. Depending on initial approximation (the dashed lines (d, e)), we obtain either confluent distribution (d) or distribution with zero-valued ordinate (e).

(25)

where ∆ h is the mean square of relative error in n(P) that is usually known from the physical experiment. With a large N one may expect a transgression of uniqueness of solution. In particular, the number of optimal sets of amplitudes is infinite for confluent distributions. Though the number of optimal solutions is finite for almost confluent distributions, the occurrence of several minima of ∆ causes a problem of physical interpretation of a solution. Probably, an appearance of the local minima of ∆ is possible even if criterion 25 is violated. However, with a small N (practically N = 1-5), the algorithm should lead to a single (global) minimum of ∆. To find a confluent distribution with positive ordinates, it is sufficient to accomplish the substitution Fi ) exp(xi). Since the algorithm copes with finding of confluent distributions, the question may be raised as to whether the barrier term in eq 22 is necessary. We deem that the use of this term is advisable at least for three reasons. First, depending on initial approximation, it allows one to find not only the confluent distributions but also those including zero-valued ordinates (see Figure 3e). Second, as numerical experiments show, it considerably widens the region of the algorithm convergence. Third, at finite values of Rs the barrier additive possesses a regularizing effect that is of special importance if N is large. It is obvious that the algorithm leads to decrease F˜i; with arbitrary N we shall have F˜i f exp(-D) as Rs f ∞ (see also the Appendix). Let us continue the consideration of numerical aspects of the algorithm. The substitution

xi ) ln(Fi)

(1 e i e N)

When equidistantly distributed on the logarithmic scale of P within the interval -1 e log10P e 1, the 11 points (see Figure 3a) have been calculated by the formula

restricts the region of the search by positive values of ordinates {Fi}. To avoid the manipulation with the constant A0, it is convenient to introduce the new searched for variables {zi} (the dimensionless, reduced energies) instead of the energies {i}

n(P) ) P[0.75/(P + exp(-1)) - 0.25/(P + exp(1))]

zi ) -ηi/RT ) -i/RT + ln A0

(24) With isotherm 9 and N ) 1 (see Figure 3c) we have MRD ) 6.067%, Fh1 ) 0.548, zj1 ) -j1/RT + ln A0 ) -1.32. The degeneracy of distribution denotes that the problem has infinite solutions. In the above example with N ) 2, the ordinates Fh1 and Fh2 can be arbitrary; it is only of importance that their sum be constant i.e., Fh1 + Fh2 ) 0.548 (see Figure 3d). The algorithm that is considered here copes with the finding of confluent distributions (see also the fourth section). While the algorithm is performed, the values of i and j for some components are nearing; thus, a sum of ordinates actually corresponds to one ordinate. The final values of ordinates depend on an initial approximation (see Figure 3d,e). Since the experimental function n(P) includes the random errors, an ideal degeneracy is hardly probable. One may say only about approximate degeneracy. We deem that the dependence of the function ∆ on N can serve as the sole measure of this degeneracy. If an addition of subsequent component with the index N + 1 does not lead to an appreciable decrease of ∆, the solution with N + 1 components can be considered as almost confluent.

In order to find the optimal values of {Fi} and {i} (or {xi} and {zi}) some iterative procedure must be used. Let us suppose that after the Ith iteration, the search produces the values {xi} and {zi}. At the next iteration (I + 1) they will receive the relative increments {βi} and {γi} and will be equal tond

x(I+1) ) (1 + βi)x(I) i i ) (1 + βi)xi

(26)

z(I+1) ) (1 + γi)z(I) i i ) (1 + γi)zi

(27)

(1 e i e N) The functions {Gi} for each value of P will also be changed

) exp[xi(1 + βi)] θ[P,(1 + γi)zi] ≈ G(I) G(I+1) i i + βiBi + γiΓi ) Gi + βiBi + γiΓi (28) (1 e i e N)

3636



problem of conditional minimization

where

∆Rs*[{βi},{γi}] w min

Bi ) (∂G(I+1) /∂βi)|βif0,γif0, Γi ) (∂G(I+1) /∂γi)|βif0,γif0 i i

N

x(I) i ,

z(I) i ,

β2i ) Sβ ∑ i)1

G(I) i

In eqs 26-28 the variables and are denoted as xi, zi, and Gi to simplify the expressions below. Let us write the partial derivatives of the functions G(I+1) in the explicit form. For the model of Langmuir 9 i and the model of random heterogeneity 11, they are identical

Bi ) xi exp(xi)θi,

Γi ) -zi exp(xi)θi(1 - θi)

where θi ) θ(P,zi) or θi ) θ(P,zi,Θ). For the patchwise model 10 the derivatives are written as

N

γ2i ) Sγ ∑ i)1 where Sβ1/2 and Sγ1/2 are the total lengths of increments. As R1 f ∞ and R2 f ∞ we have, correspondingly, Sβ f 0 and Sγ f 0. With small Sβ and Sγ, linearizing 28 gives high accuracy (I+1) ∆Rs[{Gim }] ≈ ∆Rs*[{βi},{γi}]

that is why minimization 31 will lead to the decrease of ∆Rs, i.e.

Bi ) xi exp(xi)θi Γi ) -zi exp(xi)θi(1 - θi)/[1 - (Zw/RT)θi(1 - θi)] It is convenient to consider the new function

∆*Rs[{Gim + βiBim + γiΓim}] ) ∆Rs*[{βi},{γi}] (29) where Gim ) Gi(Pm), Bim ) Bi(Pm), and Γim ) Γi(Pm). To find a minimum 15, 16 (or minimum 23 with Rs f 0), it is necessary at each iteration to search for the unconditional minimum of the function 29 with respect to the 2 × N variables {βi} and {γi}, i.e.

∆Rs*[{βi},{γi}] w min

(I+1) (I) }] < ∆Rs[{Gim }] ∆Rs[{Gim

Unfortunately, the rate of the convergence will decrease as R1 f ∞ and R2 f ∞. Thus, one needs to find some compromise by choosing the proper values of R1 and R2. This can be performed empirically by a series of trial calculations. Note that in the normal iterative process βi f 0 and γi f 0; therefore, the minima of ∆Rs,R1,R2 and ∆Rs* converge to common minimum 23. The practical calculations have shown that scaling (26 and 27) is a necessary element. When solving the problem with respect to the absolute increments, the stabilizer in eq 31 will have the form

(30)

The function ∆Rs*[{βi},{γi}] has a quadratic form; therefore at each iteration the minimization 30 will result in the unequivocal minimum of ∆Rs*[{βi},{γi}]. However, (I+1) the original function ∆Rs[{Gim }] differs from ∆Rs*[{βi},{γi}] for the latter function includes the errors of linearizing eq 28. Because of these errors, one may expect the situation when with decreasing ∆Rs*[{βi},{γi}] the increase (I+1) of ∆Rs[{Gim }] will occur; i.e., at the minimum 30 it will be

N

(I) (I+1) }] < ∆Rs[{Gim }] ∆Rs[{Gim

N

- xi, ∆zi ) z(I+1) - zi denote the absolute where ∆xi ) x(I+1) i i increments. The values of xi and zi at different i considerably differ from each other, as a consequence, the regularization effect is substantially weakened. Since the function ∆Rs,R1,R2 has a quadratic form, minimization 31 comes to solving the system of 2 × N linear equations following from the conditions

∂∆Rs,R1,R2/∂βi ) 0,

(I) ∆Rs*[{βi},{γi}] < ∆Rs*[{Gim }] )

N

N

(∆xi)2 + R2∑(∆zi)2 ) R1∑x2i β2i + R2∑z2i γ2i ∑ i)1 i)1 i)1 i)1

R1

∂∆Rs,R1,R2/∂γi ) 0

(1 e i e N) (32)

Having deduced the expressions for the partial derivatives 32, we have

This obviously means the divergence of iterations. To avoid the mentioned effect, one must resort to regularization. Its simplest variant is realized by minimizing the following function:

N

∑ i)1

N

bj,iγi + a j jβj + bh jγj ) cj ∑ i)1

aj,iβi +

(1 e j e 2 × N)

N

β2i + ∑ i)1

∆Rs,R1,R2[{βi},{γi}] ) ∆Rs*[{βi},{γi}] + R1 N

R2

∑ i)1

γ2i

w min (31)

where M

aj,i ) M-1

Θ-2(Pm)Bi(Pm)Bj(Pm) ∑ m)1 M

where R1 and R2 are the positive parameters which can be considered formally as Lagrange multipliers in the

bj,i ) M

-1

Θ-2(Pm)Γi(Pm)Bj(Pm) ∑ m)1

(33)

Analysis of Surface Heterogeneity M

cj ) M-1

∑

Langmuir, Vol. 12, No. 15, 1996 3637

or

N

Θ-1(Pm)[1 - Θ-1(Pm)

m)1

∑ i)1

Gi(Pm)]Bj(Pm) qst ) B/A

Rsxj(xj + D) a j j ) Rsx2j + R1,

b hj ) 0

where A and B vary according to models 9, 10, and 11 to give the following: for the model of Langmuir

(if 1 e j e N)

M

aj,i ) M-1

Θ-2(Pm)Bi(Pm)Γj(Pm) ∑ m)1

N

A)

M

bj,i ) M-1

Θ-2(Pm)Γi(Pm)Γj(Pm) ∑ m)1

M

cj ) M-1

(36)

Fiθi(1 - θi) ∑ i)1

N

Fiθi(1 - θi) dzi/dT ∑ i)1

B ) RT2

N

Gi(Pm)]Γj(Pm) ∑ Θ-1(Pm)[1 - Θ-1(Pm)∑ m)1 i)1

for the model of patchwise heterogeneity N

a j j ) 0,

b h j ) R2

(if N < j e 2 × N)

A)

At each iteration (I ) 0, 1, 2, ...) the system of the linear equation 33 is solved, then correction according to formulas 26 and 27 is carried out. The computational experiments show that for converging the iterations it is necessary to choose well the initial approximation. However, we find it hard to recommend any universal procedure for this purpose. The optimization using the equidistant nodes seems to be the most reasonable method for a crude estimation of {xi}. After prior optimizing of {xi}, one may proceed to the correction of both ordinates {xi} and energies {zi}. We note also that the initial xi and zi should not be equal to zero, otherwise in accordance with eqs 26 and 27 we shall have xi ) 0 and zi ) 0 for all iterations.

N

B)

Fi{[θi(1 - θi)]-1 - Zw/RT}-1(RT2 dzi/dT + ∑ i)1 θiZw)

for the model of random heterogeneity N

A)

qst ) RT2(d ln P/dT)|Θ)const

(34)

We shall suppose that the {Fi} and {zi} values are known; they are calculated from an adsorption isotherm by the above described method. To obtain qst, one must find, first of all, the derivatives of the local isotherms 9, 10, and 11 upon T: for the model of Langmuir

dθi/dT ) θi(1 - θi)(d ln P/dT - dzi/dT) for the model of patchwise heterogeneity

dθi/dT ) θi(1 - θi)(d ln P/dT - dzi/dT - θiZw/RT2)/ [1 - Zwθi(1 - θi)/RT] for the model of random heterogeneity

dθi/dT ) θi(1 - θi)(d ln P/dT - dzi/dT - ΘZw/RT2)

N

RT2 dzi/dT ) i + RT2 d ln A0/dT ) RT(ln A0 + T d ln A0/dT - zi) (37) We shall assume that A0 does not depend on i, i.e., A0(1) ≈ A0(2) ≈ ... ≈ A0(N). This approach is generally accepted most commonly (see surveys1-3 and recent publications8-10). Although vibration energy is related to the adsorptive potential, the differences i - j are usually considered to be much smaller than the average value of the potential j ≈ 1 ≈ ... ≈ N, i.e. i - j , j. Note that a low-temperature limit of A0, when a degeneracy of oscillations takes place, does not depend on adsorption energies {i} at all. Taking account of the above assumption in a case of the Langmuir local isotherm, e.g., we obtain

∑ i)1

N

qst )

N

Fiθi) )

∑ i)1

Fi dθi/dT ) 0

Fiθi(1 - θi)[RT2 dzi/dT + ΘZw] ∑ i)1

Note that to find qst according to formula 36 it is necessary to know the set {θi} for each value of pressure P. The recalculation of {zi} into the energies {i} is performed by using the relationships

The condition Θ ) const in eq 34 means that

d/dT(

(35)

N

Fiiθi(1 - θi)/∑Fiθi(1 - θi) + ∑ i)1 i)1

RT2 d ln A0/dT ) RT[ln A0 + T d ln A0/dT N

By substituting dθi/dT from above into eq 35 and by deducing the value of qst as a multiplier, one can obtain the expression

RT2 dΘ/dT ) Aqst - B ) 0

Fiθi(1 - θi) ∑ i)1

N

B)

3. Calculation of Isosteric Heat The expressions for the isosteric heat qst follow from the equation of Clausius-Clapeyron19

Fi{[θi(1 - θi)]-1 - Zw/RT}-1 ∑ i)1

∑ i)1

N

Fiθi(1 - θi)] ∑ i)1

ziFiθi(1 - θi)/

(38)

Using eqs 36 and 37, one can readily derive equations like eq 38 for isotherms 10 and 11.

3638



It is not difficult to find the asymptotic values of qst for P f 0 and for P f ∞. With the model of Langmuir we obtain

{

qr ) qst - RT d(T ln A0)/dT ) N

ηl ) -RTzl )

N

Fiηi exp(ηi/RT)/∑Fi exp(ηi/RT), ∑ i)1 i)1 Pf0 N

ηh ) -RTzh )

N

Fiηi exp(-ηi/RT)/∑Fi exp(-ηi/RT), ∑ i)1 i)1 Pf∞

Let the adsorption energies be numbered in increasing order: η1 < η2 < ... < ηN. If they differ considerably from each other, then

{

η ≈ η ) -RTz , P f 0 qr ) ηl ≈ ηN ) -RTzN, P f ∞ h 1 1 For both models 10 and 11 taking into account the lateral interactions, the limit of qr under P f 0 is the same as for the model of Langmuir; i.e., it is equal to ηl, whereas the qr limit under P f ∞ is equal to ηh + Zw. To correct the isosteric heat for the effect of multilayer adsorption, the Clausius-Clapeyron equation for pressure of saturated vapor of the liquid adsorbate is employed

L ) RT2 d ln P0/dT where L is the heat of vaporization. By using eq 14 and the condition

RT2 dΘB/dT ) {Aqst - B + [P/(P0 - P)]Θ(P)(qst L)}[P0/(P0 - P)] ) 0 we obtain

qst ) [LPΘ(P)/(P0 - P) + B]/[A + PΘ(P)/(P0 - P)]

(39)

where A and B are calculated according to the formulas above. We note that the expressions considered can be easily generalized for a case of continuous distributions because after discretization of integrals 3 and 4 all formulas will be analogous to those deduced above for discrete distributions. 4. Numerical Results and Discussion To test the algorithm, an integral isotherm was generated numerically with the known {Fi}, {i} sets (the direct problem), and then the back calculation of the {Fi}, {i} (the reverse problem) was carried out by using this integral isotherm. In Figure 4 are shown the integral isotherms and the isosteric heats obtained by solving the direct problem for three models of local equilibrium 9, 10, and 11 with the following parameters:

N ) 4, F1 ) F2 ) F3 ) F4 ) 0.25, Zw/RT ) 3.8, zi ) - 2 - 4(i - 1) (i ) 1, 2, 3, 4) Isosteric heats (see Figure 4a) are depicted on the special scale so that the calculated data would not depend on the function A0(T) as well as on the value of d ln A0/dT. One

Figure 4. Isosteric heats (a) and the integral isotherms (b) calculated for the discrete distribution shown in Figure 5a: 1, the model of Langmuir (eq 9); 2, the model of random heterogeneity (eq 11); 3, the model of patchwise heterogeneity (eq 10).

can see that both Θ and qst in the scale of ln P are the “wavy” curves, although with the local isotherms 9 and 11 the “waves” on the integral isotherms are not so distinct as with the local isotherm 10. The “waves” on the isotherms correspond to “waves” on the isosteric heats. In the models taking into account the lateral interactions, the dependencies qst(ln P) have the maxima which correspond to the adsorption sites of different kinds. In the case of the model of patchwise heterogeneity, the maxima are especially legibly being manifested. In the model of Langmuir these maxima disappear, and with increasing the ln P the reduced isosteric heat decreases monotonically from ηl/RT ≈ -z4 to ηh/RT ≈ -z1, i.e., approximately, from the value corresponding to the interaction energy of an adsorbate with the most “strong” adsorption sites to that corresponding to the adsorption energy of the most “weak” sites. If the lateral interactions take place in adsorption system, then the final value of the reduced isosteric heat is greater than ηh/RT ≈ -z1 on the magnitude Zw/RT (see the third section). It follows from Figure 4a that the dependencies qst have quite certain view for different models; therefore, the data about the isosteric heats can be very useful for qualitative characterization of adsorption systems. As to quantitative evaluation of energetic heterogeneity it, from our viewpoint, can be made more easily by using the adsorption isotherms. In Figure 5 is represented the solution of the reverse problem with the integral isotherm Θ(P) calculated according to the model of Langmuir with N ) 4 (see curve 1 in Figure 4b). There were taken 51 points distributed evenly within the interval: -18 < ln P < 1. Initial approximation was set as Fi ) 2i/[N(N + 1)], zi ) 1 - 19(i - 1)/(N - 1), i ) 1, 2, ..., N. With Rs ) 0, R1 ) 10-3, and R2 ) 10-3, the 10-20 iterations are sufficient for the function ∆[{xi}, {zi}] to reach a constant value with accuracy its representation in computer (seven decimal digits). If an isotherm is generated with maximal accuracy (seven digits), then the {Fi}, {zi} values are calculated with an accuracy of five to six decimal digits (see Figure 5a). When the isotherm includes perturbations, the accuracy decreases. For instance, in Figure 5b is shown the distribution that was


Langmuir, Vol. 12, No. 15, 1996 3639

Figure 6. The results of approximation of the isotherm of N2 adsorption on silica at 78 K6,26 (P in Torr) with N ) 2 for (1) the model of Langmuir (eq 9), (2) the model of random heterogeneity (eq 11), and (3) the model of patchwise heterogeneity (eq 10).

Figure 5. Discrete distributions calculated for the model of Langmuir (see curve 1 in Figure 4b). The distributions corresponding to unperturbed integral isotherm with N ) 4 (a, b, d) and with N ) 10 (c) are shown by the solid lines. The results obtained after perturbation of the integral isotherm with N ) 4 (b) and with N ) 10 (d) are drawn by the dashed lines. The distributions are plotted on a scale of reduced energy: z ) -η/RT.

calculated with the isotherm perturbed on 2%. For the perturbations to be introduced into the integral isotherm, it was multiplied by (1 + 0.02Y), where Y values are the quasirandom numbers distributed evenly25 on the interval [-1,1]. It is seen in Figure 5b that in spite of error in Θ(P), calculated distribution still has a good accuracy. The solution of the reverse problem looks different if N is assigned erroneously, e.g., N is equal to 10 instead of 4. Let us name each pair of the Fi, zi values a “band” by employing a spectral terminology. If perturbations are not brought into Θ(P) (the case of a degeneracy), then with N ) 10 the grouping of bands takes place (see Figure 5c). The groups of bands (Ii bands in the set with number i) are characterized by the equal energies {RTzi} and by the set of different amplitudes {Fj,i}. The sums of the Ii amplitudes, however, are equal to Fi, i.e. ∑j)1 Fj,i ) Fi. Under isotherm perturbation the splintering of the bands takes place (see Figure 5d). Most appreciably they are shifted from their accurate positions in the region of low adsorption energies. If a similar situation would arise when treating the experimental data, the interpretation of such results would be, obviously, difficult. To substantiate the choice of N in the theoretical example, it is sufficient to resort to the simplest analysis of errors. Indeed, the increase of N from 4 to 10 does almost not affect the accuracy of approximation of the integral isotherm. In the calculations with N ) 4 and N ) 10, the MRD is equal to 1.048% and 1.010%. Thus, according to criterion 25, the distribution shown in Figure 5d is almost confluent. Without an appreciable loss of accuracy one may decrease N from 10 to 4 to obtain the result shown in Figure 5b. However, when treating real experimental data the correct choice of N is possible only in the situation (25) Forsithe, G. F.; Malcolm, M. A.; Moler, C. B. Computer Methods for Mathematical Computations; Prentice-Hall; Englewood Cliffs, NJ, 1977.

if it is known a priori that the searched for distribution is a discrete one and a model of local equilibrium used is accurate enough. Unfortunately, such information as a rule is open to question. Let us address, e.g., the adsorption isotherm of nitrogen on maximally hydroxylated silica26 (see Figure 6) that was treated earlier6 by using the algorithm IRA destined for evaluation of continuous distributions. Since, the continuous distributions6 represent well-resolved peaks, one can expect the results obtained on the basis of alternative hypotheses, which imply either continuous or discrete distribution, should be comparable. In Figure 7 calculated continuous and discrete (N ) 2) distributions are represented for each of three models with Zw/RT ) 3.8. The computations were carried out without correction for the effect of multilayer adsorption. To estimate Θ(P) in the right-hand part of eq 11, the na ) 13.13 µmol/m2 calculated with the Langmuir isotherm was used (see below). As seen, the positions and the ratios of heights of peaks corresponding to the similar models are in good agreement with each other, but only in calculations with the Langmuir local isotherm does the approximate integral isotherm provide for satisfactory approximation (MRD ) 1.825%) of experimental data (see Figure 6). The parameters naF1 ) 8.23 µmol/m2, naF2 ) 4.90 µmol/m2, z1 ) 2.561, and z2 ) -3.070 have resulted from optimization with the model of Langmuir. The calculated adsorption capacity na ) 13.13 µmol/m2 in this case is close to that obtained6 for continuous distribution: na ) 15.5 µmol/m2. At the same time models 10 (MRD ) 15.032%) and 11 (MRD ) 4.845%) evidently are not adequate for the experiment. However, one needs only to increase N and any of the models will give a satisfactory approximation. In Figures 8-11 are shown the integral isotherm (see Figure 8), the histograms, and the integral energetic distributions (see Figures 9-11) calculated with N ) 6. When calculating the discrete distributions, the continuous distributions earlier6 obtained have been used for setting the initial approximation. All discrete distributions are almost confluent, therefore depending on the initial approximation the algorithm leads to different minima of ∆, but as the computational experiments show, the distributions depicted in Figures 9-11 correspond to the most deep minima of the objective function ∆. With N ) 6 only the model of patchwise heterogeneity 10 gives appreciable error (see Figure 8) for the approximation (MRD ) 3.054%), whereas the local isotherms 9 and 11 give MRD ) 0.128%, na ) 14.17 µmol/m2 and MRD ) 0.057%, na ) 13.59 µmol/m2. Note that the error 0.057% has not been reached earlier both with N ) 2026 and with (26) House, W. A. J. Colloid Interface Sci. 1978, 67, 166.

3640



Figure 9. Histogram (a) and the integral distribution function (b) calculated for the Langmuir model (eq 9) with N ) 6. The smooth curve is a continuous distribution calculated earlier6 (see Figure 7a).

Figure 7. Continuous6 and discrete distributions calculated from the integral isotherm shown in Figure 6 for (a) the model of Langmuir (eq 9), (b) the model of patchwise heterogeneity (eq 10), and (c) the model of random heterogeneity (eq 11).

Figure 8. Same as in Figure 6 with N ) 6.

N ) 61.6 Thus, the error of approximation is not a sufficient criterion to decide the dilemma between discrete and continuous distributions. We note also that it is impossible to estimate the constant w as a posteriori result by using an integral isotherm until we have no strict information about N. Indeed, with a large N an integral isotherm can be approximated both by using the Langmuir local isotherm (w ) 0) and by using the isotherms 10 and 11 with Zw/RT ) 3.8. When normalizing the histograms plotted on the basis of formulas 7 and 8 one needs to know the value of za. The calculations with the Langmuir isotherm lead to F1 = 0 (see Figure 9); therefore one may assume za ) z1. However, the algorithm generally does not allow one to find za. To normalize the histograms in the case of isotherms 10 and 11, we have chosen za almost arbitrarily by using the formula za ) 2z1 - z2 (see Figures 10 and 11). Note that inclusion of any additional point with zero ordinate Fi ) 0 does not influence MRD.

Figure 10. Same as in Figure 9 with the model of patchwise heterogeneity.

As seen, the histograms calculated qualitatively agree with the continuous distributions. The marked difference consists of the fact that the continuous distributions have unconditional maxima in the region of low adsorption energies. However, at least for the isotherms 10 and 11 these maxima, in our opinion, can arise from the use6 of zero-order stabilizer

Ω[F] )

∫zz F2(z) dz b

a

If an integral isotherm is weakly sensitive to the values of F(z) at the ends za and zb, the algorithm prefers to decrease F(za) and F(zb). In any case the histograms calculated do not contradict the conclusions made earlier6 about occurrence of the sites


Langmuir, Vol. 12, No. 15, 1996 3641

Figure 12. Isosteric heat calculated for three models of local equilibrium with N ) 2. The numeration of curves corresponds to Figure 6. Curves 1′-3′ are calculated with correction for effect of multilayer adsorption (see eq 39).

Figure 11. Same as in Figure 9 with the model of random heterogeneity.

of two different kinds on the surface. Though, to ascertain details of the energetic distribution more accurately, it is desirable to have an integral isotherm measured at temperature lower than 78 K. To reveal the effects related to the lateral interactions, it is of interest to consider the behavior of isosteric heat. In Figure 12 is shown the qst(n) dependence that has been calculated by Aristov and Kiselev27 from the series of the adsorption isosters of nitrogen on maximally hydroxylated silica. It strictly corresponds to the system under consideration. To plot qst(n) by using any model from those above considered, one needs to find, according to eq 38, the only constant

ξ ) RT ln A0 + RT2 d ln A0/dT The value of ξ may be determined by operation of parallel shift providing for the best fit of theoretical and experimental dependencies qst(n). When N ) 2 this operation in case of the Langmuir local isotherm leads to ξ ) 12.5 kJ/mol. The same value was used for calculation with isotherms 9, 10, and 11, both with N ) 2 and with N ) 6. Unfortunately, the frequencies νx, νy, and νz in eq 13 are unknown. Thus, even though the constant ξ is available, one cannot find the constant A0. However, when using eqs 12 and 13, one can make the lower and upper estimations for A0. The low-temperature limit (T f 0, d ln A0/dT f (5/2)T-1) and the high-temperature limit (T f ∞, d ln A0/dT f - (1/2)T-1) give

RT ln A0 )

{

ξ - (5/2)RT, T f 0 ξ + (1/2)RT, T f ∞

From here we have 10.9 kJ/mol < RT ln A0 < 12.8 kJ/ mol. The magnitude RT ln A0 ) 10.5 kJ/mol (A0 ) 1.122 × 107 Torr) obtained by Jaroniec28 is very close to the lower boundary of this interval, i.e., to the low-temperature limit. (27) Aristov, B. G.; Kiselev, A. V. Zh. Fiz. Khim. 1964, 38, 1984. (28) Jaroniec, M. Surf. Sci. 1975, 50, 553.

Figure 13. Same as in Figure 12 with N ) 6.

In the absence of any information about νx, νy, and νz for estimation of A0, it is reasonable to accept d ln A0/dT ) 0. Having calculated A0 by assuming d ln A0/dT ) 0, one can, in accordance with eq 37, recalculate {zi} into {i} and plot the distributions on the scale of the dimensional energy (see Figures 7 and 9-11). In Figure 12 is shown the dependence qst(n) calculated with N ) 2. It follows from Figure 12 that the lateral interactions lead to an appearance of maxima on qst(n). Hence, in this case the Langmuir local isotherm yields the best result not only for the integral isotherm but also for the isosteric heat. The increase of N from 2 to 6 leads to better agreement of the theoretical and experimental qst(n) (see Figure 13). Note that the dependencies qst(n) calculated with isotherms 9 and 11 almost coincide. From here one may conclude that different models giving close results when approximating integral isotherms give also close results when calculating the isosteric heat. Model 10 approximates both the integral isotherm and qst(n) insufficiently, although one can expect that any model 9, 10, or 11 must give one and the same result if N f ∞. Both models 9 and 11, in our opinion, satisfactorily describe qst(n), with an exception of the region of small n. However, one must recognize that when calculating qst from a series of isotherms, appreciable errors may arise. They can be particularly large at small Θ(P). For instance, to calculate qst at Θ = 0.02-0.05 Aristov and Kiselev27 have used only two isotherms. In this case the error in qst, from our viewpoint, can attain =4-5 kJ/mol. Apparently, the direct

3642


calorimetric measurements would give more accurate and reliable information about isosteric heat.


Appendix 1. When Rs f 0, the minima of ∆ and ∆Rs coincide. Let us denote N

5. Conclusions When treating the experimental data by using considered algorithm, it is most important to set properly the number N. The sole criterion for choice of N is experimental error. If N is not large enough to decrease ∆ to the value of mean square of error in a measured isotherm, one can try to increase N. The small value of ∆, of course, does not prove the model adequacy. Moreover, if experimental data are related to high temperatures, then the magnitude of error in approximating an integral isotherm does not allow one to decide the dilemma between discrete and continuous distributions. To ascertain the nature of energetic distribution it is necessary to have thermodynamic data obtained at sufficiently low temperatures. Under small temperatures in the case of discrete distribution both an integral isotherm and a dependence of isosteric heat on pressure have typical steps (or “waves”) the number of which is equal to N. If an integral isotherm is a smooth function, one cannot find the energy of the lateral interactions w as an a posteriori result of treating the isotherm, for different models irrespective of w give close results of approximation as N f ∞. However, when an integral isotherm represents a stepwise function (low temperatures), the constant w, probably, can be estimated by using the simplest optimization principles. In any case, besides the optimization techniques it is desirable to have additional physical data about the system under consideration such as, e.g., crystallographic or spectral information concerning the surface. Furthermore, to verify a calculated distribution, other thermodynamic characteristics of the adsorbed phase such as integral and differential heats, heat capacities, and entropy should be described in addition to an integral isotherm.

V[{Fi}] )

[ln(Fi) + D]2 ∑ i)1

If {Fhi}, {zji} correspond to a minimum of ∆, the insertion of any other sets {Fi}, {zi} in ∆ can only increase the ∆. In particular, if {F˜i}, {z˜ i} are the sets corresponding to a minimum of ∆Rs and λ is a small positive constant, one may write the following inequalities:

∆[{Fhi},{zji}] < ∆[{Fhi},{zji}] + RsV[{F˜i}] < ∆Rs[{F˜i},{z˜ i}] e ∆Rs[{Fhi+λ},{zji}] ) ∆[{Fhi+λ},{zji}] + RsV[{Fhi + λ}] Let us choose so small a λ that [ln(λ) + D]2 > [ln(λ + Fhi) + D]2 (1 e i e N), then

∆[{Fhi},{zji}] < ∆Rs[{F˜i},{z˜ i}] < ∆[{Fhi+λ},{zji}] + RsV(λ) Let us choose an arbitrary small positive constant γ. Since ∆ is a continuous function, one can find such λ that ∆[{Fhi+λ},{zji}] < ∆[{Fhi},{zji}] + 0.5γ. We always have a possibility of taking Rs from the relation Rs ) 0.5γ/V(λ) to obtain

∆[{Fhi},{zji}] < ∆Rs[{F˜i},{z˜ i}] < ∆[{Fhi},{zji}] + γ 2. At the minimum point of ∆Rs the values of ∆Rs and ∆ are the monotonically increasing functions of Rs while V is a monotonically decreasing function of Rs. Let ∆Ri s, ∆i, and Vi be the values of above functions at the minimum point corresponding to Rs ) Ris. If R2s > R1s then

∆1 + R2s V1 > ∆R2 s ) ∆2 + R2s V2 > ∆2 + R1s V2 > ∆1 + R1s V1 ) ∆1Rs Hence, ∆R2 s > ∆R1 s, R2s (V1 - V2) > ∆2 - ∆1 > R1s (V1 - V2), - R1s )V1 > (R2s - R1s )V2, V1 > V2, ∆2 > ∆1.

(R2s

LA950181W

Numerical Method for Analysis of Surface ... - ACS Publications

Recommend Documents