Abstract
Transcendental functions are a fundamental building block of science and engineering. Among them, a relatively new function denominated as Lambert is highlighted. The importance of such function relies on the fact that it can perform novel isolation of variables. In this work, we propose two accurate piece-wise approximate solutions, one for the lower branch and another one for the upper branch, respectively. The proposed analytic approximations are obtained by using the power series extender method (PSEM) in combination with asymptotic solutions. In addition, we will compare some published approximations with our proposal, highlighting our advantages in terms of significant digits and speed of evaluation. Furthermore, the approximations are validated by the successful simulation of a problem of economy and other acoustic waves of nonlinear ions.
1. Introduction
Transcendental equations are those containing transcendental functions with the variable to be solved for, which have the particularity that can not be expressed in terms of a sequence of algebraic operations [1–4]. This is the case of the Lambert Function [5, 6]. The Lambert W Function as today known, denoted by , is a multivalued function that can be defined as the inverse of the following transcendental function [5–9]:resulting inwhere can be any real or complex number, and is known as Lambert Function.
The Lambert W Function () is a multivalued function (also called multibranched function); branches considering only real numbers of are plotted in Figure 1. The two real branches, and , can be identified. is called the upper branch (also named the principal branch) and satisfies the condition . This branch can be split into two branches, and . is called the lower branch and satisfies the condition . Both branches are jointed at (, ). Then, on the interval , there are two possible values for , the one for and the second for .
The origin of (2) begins in 1758 when Johann Heinrich Lambert solved the next trinomial equation [5, 6] by giving a series for in powers of .
Some years later, in 1783 Leonhard Euler presented (3) into a symmetrical form as follows [5, 6]:
which is developed in series expansion as
and can be expressed as
This series is called the Tree Function and is related to the Lambert W Function (2) because [5, 6]. The Lambert W Function is also called the Omega function or Product Log function [5, 13]. However, the inverse of (1) was first described by Polya and Szego in 1925 [14]. Later, during 1949 up to 1959, Edward Maitland Wright performs significant contribution to the solution of the transcendental equation (1) by means of his effort to solve Delay Differential Equations (DDE) [15, 16]. At the same time as Wright, in 1950, N. D. Hayes [17] collaborated in the solution of (1) with his work on transcendental equations. In 1973, Crowley et al. [18] presented an algorithm (called the 443-algorithm) for the solution of that transcendental equation; an updated version of the 443-algorithm (called the 743-algorithm) was given in 1995 [19]. Different works have been published from 1758 along the years regarding the solution of (1); however, it was until the work of Corless et al. published in 1996 [5] that is given to be known worldwide. Corless et al. [5] made a revision of related works and potential applications and provide approximate representations and numeric algorithms of evaluation to the function. At the same time, this function was included in Maple CAS software [5, 20].
Since then, many works have appeared with reviews and significant contributions to the approximation of the oriented to diverse applications; some relevant works are presented by Corless et al. [21, 22], Brito et al. [6], Abdolhossein et al. [23], Fukushima [24], Marko [25], Fredrik [9], Jeffrey et al. [26], Jeffrey [27], Kheyfits [28], Darko [7], Dence [29], and Mezo et al. [30]. Moreover, a conference was created and celebrated in honour of the Lambert W Function in 2016. Lambert is today used in all fields of science where transcendental functions as (1) are implicit in equations to solve. In general, can be applied in Physics, Mathematics, Chemistry, Engineering, and so on [5, 6, 8]. Some particular applications are biochemical kinetics [25, 31], electromagnetics (capacitor, transmission lines, etc.) [32], semiconductor (diode) circuits [33, 34], photocells [35], acoustics [36], optics [37], particle physics [38], general relativity [39], geophysics [40], risk theory [41], technological systems [42], information theory [43], combinatorial [5, 44], iterative exponentiation [5, 45, 46], graph theory [47], oscillation theory [13], time-delay dynamics [48–51], pedagogical [5], physiology [52, 53], hydrology [54, 55], hydraulic [56], materials and transport research [57–59], electrochemistry [60], microfluidics [61], statistical mechanics [62], etc.
The formal solution of the Lambert W function can be expressed as [56]
This expression clearly requires an infinite iterative calculus, which is not feasible to compute.
From the numerical and analytical point of view, there are several methods to compute transcendental functions like this; for example, Chapeau-Blondeau [10] and Corless et al. [5] proposed to evaluate the Lambert W branches using approximate analytical expressions and increasing accuracy with the Halley numerical method [35, 63–65]. Other proposals to obtain an analytic approximation of the Lambert W function are [11, 66–68]. Fukushima proposed an algorithm [24] to obtain the numerical solution for and by means of an expansion around the zero branch point. Moreover, Chapeau-Blondeau [10], de Brujin [67], and Karamata [68] proposed series expansions to evaluate branch .
In this work, we propose to apply power series extender method (PSEM) [69] to calculate approximate solutions for both branch Lambert . PSEM is a novel tool that is able to produce easy computable, highly accurate approximations without requiring to perform: integrals, solution of differential equations or inclusive without requiring the existence of a perturbative parameter. What is more, the approximate solutions exhibit a large domain of convergence.
Section 2 presents the basic procedure to apply PSEM. The procedure to obtain the approximation of lower and upper branch of Lambert is presented in Section 3. Next, two case studies are presented in Section 5. Later, Section 6 presents a discussion about the main results of this work. Finally, the concluding remarks are presented in Section 7.
2. Basic Procedure of PSEM
In broad sense, a nonlinear differential equation can be expressed as follows:
having as boundary condition,
where and are a linear operator and a nonlinear operator, respectively, is a known analytic function, is a boundary operator, is the boundary of domain , and denotes differentiation along the normal drawn outwards from [70].
Next, we express the solution of (8) as a Taylor expansion:where are the coefficients of the power series and is the expansion point.
Now, in [69, 71], it is proposed that the solution for (8) can be written as a finite sum of functions asor
where are constants to be determined by PSEM, are arbitrary trial functions, and and are the orders of approximations (11) and (12), respectively. We will denominate (11) and (12) as a trial function (TF). Next, we calculate the Taylor expansion of (11) or (12), resulting in the power series:
respectively, where Taylor coefficients are expressed in terms of parameters . Finally, we equate/match the coefficients of power series (13) or (14) with the correspondent ones of (10) to obtain the values of and substitute them into (11) or (12) to obtain the PSEM approximation [69].
It is possible to enrich PSEM by means of introducing cancellation points (CP). For every cancellation point several equations depending on the number of desire derivatives to satisfy are constructed. Such new equations are included in PSEM process to complement the system of equations in order to obtain the coefficients (). The cancellation equations are constructed by evaluating the TF and its derivatives at some strategically selected CP to increase the accuracy within a region. It is important to note that every introduced cancellation point and its derivatives would replace the corresponding superior order equations in the regular PSEM procedure.
The TF can be composed by combining different type of functions as long as the Taylor series exist. Note that regrouping terms of the -powers are valid for (13) and (14), because functions were chosen analytic in a well-defined domain in the independent variable , whereby, the correspondent Taylor series are convergent for such values of [72, 73]. By restricting the values of (8) to the mentioned domain of convergence, the sum of the Taylor series of is also convergent [72, 73]. In the same fashion, the quotient of two analytic functions series (see (12)) in is also analytic, whenever the denominator is different from zero at . It is important to notice that PSEM convergence greatly depends on the proper selection of the trial function. Then, it is necessary that the proposed TF can potentially describe the qualitative behaviour of the solution of the nonlinear problem.
3. PSEM Approximation Procedure for Lambert Function
In this section, we will present the approximations for both branches of Lambert function by using the PSEM procedure. Each branch will be divided into some regions in order to reduce the complexity of the approximations increasing in consequence with the speed of evaluation. In addition, all the numerical PSEM approximations of this work are obtained considering at least four significant digits. Readers can find in Appendix A the proposed numerical approximations. It is important to remark that all the power series expansions of this work were obtained using special transformations of and applying the built-in package “MultiSeries” of Maple 18.
3.1. Lower Branch Treatment
The lower branch will be segmented into three pieces. For the first approximation, we propose
where constants and will be determined by the PSEM procedure. In order to obtain an accurate approximation, we propose to obtain build 4 equations for , evaluating them at and other 5 equations for , evaluating them at , and equating with the same derivatives of (15). Then, by solving the resulting equation system for the unknown constants we obtain the approximate solution (A.1) (see Appendix A). The second and third approximations are expressed asandwhere
It is important to note that is obtained by calculating the order two power series expansion at forand transforming (19) using , we obtain (18).
In summary, the lower branch is approximated (see (A.1)-(A.3) in Appendix A) by
3.2. Upper Branch Treatment
The upper branch is divided into four segments. First, we propose the following transform:Then, the proposed trial function to apply PSEM to (21) iswithwhere is obtained by calculating the order four power series expansion at for (21). Then, the first approximation for the upper branch of Lambert is
Next, the TF for the next two approximations areand
For the last approximation, we propose the following transform:
Then, we calculate the order 5 power series expansion of (27) and replace ; it results in
In summary, the upper branch is approximated by
Numerical PSEM approximations , , and are (A.4), (A.5), and (A.6), respectively (see Appendix A).
4. Iterative Scheme to Improve Accuracy
In [74] an easy computable iterative algorithm was proposed to increase the accuracy of approximation. The iterative formula is given bywhere is number of iterations, isand
The error term for this iterative scheme is [74, 75]. In general terms, lower and upper branch were approximated by PSEM in order to obtain more than 4 significant digits. Then, after a single pass of 30 we obtain more than 16 digits as depicted in Figure 2. In the same fashion, Figures 3 and 4 show more than 15 digits for the upper branch approximations. Therefore, our scheme is compatible with the double float (64-bit IEEE 754) format of C/C++, which provides 15 significant digits. Furthermore, the same figures presented the significant digits of other works from literature after a singles pass of (30), causing our proposal to be the only one that keeps more than 15 digits for the complete domain of Lambert . If more significant digits are required, user can add more iterations to the scheme. The significant digits were calculated usingThe number of employed digits extended to 50 (in Maple CAS software) during the simulation in order to avoid rounding errors during the evaluation of approximations. It is important to notice from Figure 4 that our approximation presents the higher number (more than 23) of significant digits for than all the other approximations.
4.1. Computation Time Measurement
A computation time (CPU time) comparison among our approximation and other ones from literature was performed using Fortran 77/90, using as compiler the gcc 5.4.0 with level 3 of optimization. The computer used for the simulations was an Intel Core Pentium(R) CPU 2117U × 2 running at 1.80GHz under Linux Ubuntu 16.04. In order to circumvent the operative system instabilities, the time measurement was determined by the average time of 10 million of evaluations for each point. It is important to note in Figure 5 that our approximation presents, in general, the best CPU time behaviour in comparison with other reported approximations: [7, 10] for lower branch and [7, 11, 12] for the upper branch. However, for some regions our approximation exhibits a similar or slightly better CPU time. Further work is necessary to obtain faster approximations compatible with double float standard as the one reported [24], which presents a precise and fast computation numerical algorithm to evaluate that exhibits 50 nsec of average time. The key advantage of [24] to obtain such speed is that authors do not use transcendental functions to evaluate Lambert . However, our approximation is simple piece-wise formulas that are straightforward to implement in Fortran (see Appendix A), C/C++, Java, Basic, among many others.
(a) C.T.
(b) C.T. for
5. Case Studies
In this section we present two applications of the Lambert function . The first case of application is to economy and second one is an interesting application to plasma physics.
5.1. Economic Order Quantity with Perishable Inventory
In [76] a problem of Economic Order Quantity (EOQ) with perishable inventory related to the supply chain design of melons was reported. After some assumptions the authors of [76] obtain the optimal transfer batch quantity aswhere the batch transfer cost is [$], the batch transfer time is [h], the time in the cold chain is , the value of the melons at picking is , the deterioration rate at a field temperature at 30°C is , the deterioration rate in the cold chain is , and the picking rate is .
As an example obtained from [76], we assume the values for $7, h, h, $75, days, per day. Allowing the picking rate to vary between 1 and 120 cartons per hour producing an optimal transfer batch quantity, , which is shown in Figure 6, where we have plotted (34) using our proposal and the exact numerical Lambert based solution from Maple 18.
5.2. Application on Extremes of the Sagdeev Pseudo-Potential
After several assumptions and manipulations, the electrostatic potential solution [77] is obtained asand
where the first and second roots correspond, respectively, to the maximum and minimum in the pseudo-potential, in terms of both real branches of Lambert . and are constants of the phenomenon.
The exact ( and ) and approximate ( and ) contour plot for different values of is presented in Figure 7(b). In addition, to find the locus of the conjugation points (see Figure 7(a)), we must insert in (35) and (36).
(a) Contour plot of the roots for and
(b) Locus of conjugation points
The conjugation points separate the (upper) region of the root , corresponding to the maximum in the pseudo-potential, from the (lower) region of the root , corresponding to the minimum in the pseudo-potential.
6. Discussion
PSEM procedure produces compact and high accurate approximations as depicted in Figures 2, 3, and 4. Furthermore, as a reference for convergence of PSEM, we use the root mean square error defined aswhere means relative error and and represent the valid region of the approximation. The integration procedure will be performed using Simpson’s rule 1/3 taking 100 samples for each valid interval of the approximations (20) and (29). Additionally, in order to calculate (37), we consider the exact value of Lambert function as the numerical values obtained from the built-in command from Maple 18.
Figure 8 depicts the error for each branch approximations, using several -th trial functions (TF) increasing gradually the number of constants accordingly with Tables 1 and 2. Such figure shows how as the number of constants increments the error reduces its value significantly showing the good convergence of the PSEM procedure.
As it is expected, the cases study shows how our proposal can be applied to solve practical problems in sciences and engineering. In fact, the simplicity of the piece-wise approximations allows easily implementing the approximations within computer languages as Fortran, C++, among many others. It is also worth mentioning that the computation time of our proposal was compared against other reported approximations from literature presenting the advantage of our proposal. Further work is required to reduce even more the computation time in order to compete with numerical algorithms as the one reported in [24].
7. Concluding Remarks
In this article piece-wise approximations for both branches of Lambert were proposed. The approximations were obtained using PSEM in combination with asymptotic approximations, resulting in handy and highly accurate expressions in comparison with other ones from literature. In fact, we apply the approximations to practical problems in the fields of economic and physics, exhibiting highly accurate results. It is important to highlight that the coupling of PSEM with asymptotic solutions was successful to improve the accuracy in the vicinity of the expansion point and keeping the high accuracy of the original asymptotic approximation for large values of the domain. Additionally, we prove that it is possible to improve the accuracy using an easy computable iterative scheme based on asymptotic methods. The computation time was measured resulting that our proposal is fast in comparison with other approximations from literature. Even more, the significant digits from our proposal are fully compatible with the double float format from C/C++ concluding that our approximations can be applied to practical problems as the study cases of this work.
Appendix
A. Numerical PSEM Approximations
The numerical approximations obtained using PSEM with parameters from Table 3 are for the lower branch:and
Next, the approximations for the upper branch areand
B. Fortran Code to Compute the CPU Time
!===================================================================
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Routine for measure CPU time of the PSEMapproximations for Lambert W
real:: start, finish, tt
integer jend,j
integer8 i,ii
real8 z,x,w,wmin,wmax,dw
real8 V,A, Q,F
real8 k0,k1,k2
real8 k3
ii=10000000
jend=120
wmin=-15.0d0
wmax=15.0d0
dw=(wmax-wmin)/dble(jend)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
write(,"(a20,a20,a25)") "z","W","Tp nS"
do j=0,jend
x=wmin+dwdble(j)
w=x
call cpu_time(start)
x=wexp(w)
do i=0,ii
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Lower branch LambertW!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!! W-1,3(x)
if(-0.1<=x.and.x<0.and. w<-1) then
k0 = (1.01999365162218+(-12.6917365519443-45.1506015092455x)x)x/(1+&
(-22.9809693297808+(-104.692066099727-95.2085341727207x)x)x)
k1 = log(-x)
k2 = k1-log(-k1)+log(-k1)/k1
V = k0+k2
end if
!!!!!!!!!!!!!!!! W-1,2(x)
if(-0.34<=x.and. x<-0.1.and. w<-1) then
V = (-1362.78381643109+(-1386.04132570149+(11892.1649836015+
16904.0507511421x)x)x)x/(1+ &
(251.440197724561+(-1264.99554712435+
(-5687.63429510978-2639.24130979048x)x)x)x)
end if
!
!!!!!!!!!!!!!!!! W-1,1(x)
if(-0.3678794411714423<=x.and.x<-0.34.and. w<-1) then
V= (-7.874564067684664+(-63.11879948166995+(-168.6110850408981-&
150.1089086912451x)x)x)x(x+0.3678794411714423)/(1.0+&
(15.97679839497612+(98.26612857148953+(293.9558944644677+&
(430.4471947824411+247.8576700279611x)x)x)x)x)-1.0
end if
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Upper branch!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!!!!!!!!W0,1(x)
if(-0.3678794411714423<=x.and.x<1.and. w>=-1) then
k1 = sqrt(1.+2.7182818284591x)
k2 = 0.33333333333333+.70710678118655/k1-0.058925565098880k1+&
(x+0.36787944117144)(0.050248489761611+(0.11138904851051+
0.040744556245195x)x)/(1.+ &
(2.7090878606183+(1.5510922597820+0.095477712183841x)x)x)
V = -(k2-1)/k2
end if
!
!!!!!!!!!!!!!!!!W0,2(x)
if(1<=x.and.x<40) then
k1 = 1+(5.950065500550155+(13.96586471370701+(10.52192021050505+
(3.065294254265870+0.1204576876518760x)x)x)x)x
V = 0.1600049638651493log (k1)
end if
!
!!!!!!!!!!!!!!!!W0,3(x)
if(40<=x.and.x<20000) then
k1 = 1+(-3.16866642511229e11+(3.420439800038598e10+ &
(-1.501433652432257e9+(3.44887729947585e7+(-4.453783741137856e5+ &
(3257.926478908996+(-10.82545259305382+(0.6898058947898353e-1+ &
0.4703653406071575e-4x)x)x)x)x)x)x)x)x
V = 0.9898045358731312e-1log(k1)
end if
!
!!!!!!!!!!!!!!!!W0,4(x)
if(20000<=x) then
k1 = 1/(1+log(1+x))
k2 = 1/k1
k3 = log(k2)
V = k2-1-k3+(1+k3+(-1/2+(1/2)k32+(-1/6+(-1+(-1/2+
(1/3)k3)k3)k3)k1)k1)k1
end if
!
! Fritsch's formula
A = log(x/V)-V
Q = (2(1+V))(1+V+(2/3)A)
F = A(Q-A)/((1+V)(Q-2A))
z = V(1+F)
continue
continue
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
enddo
call cpu_time(finish)
tt=((finish-start)/10000000)1000000000
write(,"(0pf25.15,0pf25.15,0PF15.5)")x,z,tt
enddo
end program lambertw_tiempos
Data Availability
All data generated or analysed during this study are included within this research article.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
Authors wish to thank, M.TN. César Colocia García, Eng. Francisco Martínez Barrios, Roberto Ruiz Gomez and Gerardo César Vélez López for their support. The second author is currently a postdoctoral fellow researcher at National Institute for Astrophysics, Optics and Electronics and expresses his gratitude for the support received. He also acknowledges the financial support from the Secretary of Public Education of México (SEP) through Grant 2703 E476300.0275652.