The Lambert function, the quintic equation and the proactive discovery of the Implicit Function Theorem

One of the problems on which a great deal of focus is being placed today, is how to teach Calculus in the presence of the massive diffusion of Computer Algebra tools and online resources among students. The essence of the problem lies in the fact that, during the problem solving activities, almost all undergraduates can be exposed to certain “new” functions, not typically treated at their level. This, without being prepared to handle them or, in some cases, even knowing the meaning of the answer provided by the computer system used. One of these functions is Lambert’s W function, undoubtedly due to the elementary nature of its definition. In this article we introduce W, in a way that is easy to grasp for first year undergraduate students and we provide some general results concerning polynomial-exponential and polynomial-logarithmic equations. Among the many possible examples of its applications, we will see how W comes into play in epidemiology in the SIR model. In the second part, using more advanced concepts, we motivate the importance of the Implicit Function Theorem, using it to obtain the power series expansion of the Lambert function around the origin. Based on this approach, we therefore also provide a way to obtain the power series expansion of the inverse of a given smooth function f (y), when it is assumed that f (0) = 0, f ′(0) 6= 0, aided by the computational power of Mathematica®. Basically, in this way, we present an alternative approach to the Lagrange Bürman Inversion Theorem, although in a particular but relevant case, since the general approach is not at an undergraduate level. A number of good references are [1, pp. 23–28] and [2], where the Lambert function is applied. Finally, these skills are used to take into consideration the particular quintic equation in the unknown y presented by F. Beukers [3]. Namely, we consider x(1 + y)5 − y = 0 as an example of an equation for which the power series representation of one of its real solutions is known, calculating, with the same method used for the Lambert function, the first terms of its power series representation.


Introduction
T he Lambert function, whose introduction historically dates back to 1758 Johann Heinrich Lambert's contributions [4] and later studied by Leonhard Euler [5] and Francesco Carlini [6], was, in some sense, rediscovered at the beginning of the second half of the 90s in a number of famous articles [7][8][9]. Since then, researchers' attention has been renewed, taking into account the use of this function in a wide range of mathematical modelling applications. As a matter of fact, searching for the "Lambert function" in the MathSciNet database returns 141 entries, while in arXiv we count 59 entries (as of February 2021), testimony to the great appeal that this function has and the multiplicity of uses and applications across all scientific fields. For this reason, in this article we do not attempt to reference all of them, quoting, for instance, [10], which, in turn, refers to four contributions that indicate the physical applications (entries [2,[11][12][13] therein) and [14], where thirty application areas are listed. Taking into account the fascination generated by the cosmological topic, in this paper we prefer to limit ourselves to quoting the [15] article, where the Lambert function is used to describe the Event horizon of a Schwarzschild black hole, the recent interesting didactic note on the exponential function [16] and, in view of the impact caused by the worldwide spread of the Covid-19 virus, when calculating the terminal effects of an epidemic, looking at Remark 4.
Our approach places emphasis on the symbolic and operational aspects related to the Lambert function taking advantage of Mathematica ® (of course, any other Computer Algebra system can be used), presenting a didactic path that can be adopted at the end of the first year, or at the beginning of the second year, in an Advanced Calculus course, leading the student to discover the impact of the Implicit Function Theorem in the determination of the inverse of a given function. In our opinion, the ease of access to Computer Algebra resources means that any instructor must be adequately updated in this regard, in order to fully assist the students in their learning process. We firmly support the claim that the Lambert function is an elementary function and, for this reason, it should be included in Calculus courses, since this function often occurs in common problem solving, following what was first stated by [7] and later by [11,13,17]. We do very much hope that our contributions may help advance this pursuit. It is also worth mentioning that the Lambert function is categorised among elementary functions on the Wolfram function website. As a matter of fact, as pointed out in [18], only elementary notions are used and needed to introduce the Lambert function. Indeed, it is formally the same as defining inverse trigonometric functions, where we take the inverse function of a given piecewise monotonic smooth function f : I → f (I) = J over the I, interval, which identifies the new function f −1 : J → I = f −1 (J). These notions are commonly covered in a first year undergraduate course, or possibly at the College level. In this regard, we see a close analogy with the use of the Euler Gamma function, in the sense that the definition of these functions depends on basic tools, while Special functions are usually introduced relying on Differential Equations.
The great applicability and benefits of the Lambert function are attributed to the fact that it naturally appears when solving a wide range of transcendental equations involving exponential and logarithmic functions. As a matter of fact, the process that leads to the definition of the Lambert function is connected to the search of roots of the simplest exponential-polynomial equation.
Lastly we highlight two facts: refencing Wikipedia or the Wolfram Mathematical function website , as well as many other internet resources, any student can obtain a great deal of information on the Lambert function, in spite of the fact that the latter is very seldom introduced by instructors and, moreover, not often referenced in the majority of introductory textbooks. What's more, the Lambert function is implemented in all Computer Algebra systems, as per trigonometric, hyperbolic, exponential or logarithmic functions, normally handled in Calculus courses.
We are in agreement with the emphasis that was already given by [17] towards experimental mathematics, while at the same time trying to outline a path that leads the student towards the conscious use of Computer Algebra. To further emphasise what we mean, let's look at a simple illustrative example, comparing the equations For a student, there is nothing strange when considering the Equation (1a), while the Equation (1b), instead, provides no evidence about its possible symbolic treatment. It is generally accepted that the fundamental (i.e., in the first quadrant) solution to (1a) is x = arcsin 1 3 and the approach to describe the whole family of solution of (1a) is a standard matter. Maybe later, depending on the level of the course, it is shown that the fundamental solution is the sum of a (beautiful) infinite series: Since, when |y| < 1 it is possible to show, using the McLaurin series, that 2n n The (2) and (3) formulas are obtained using the derivatives of the proper local inverse of the sine function. Our aim is to show that the same approach, which can be followed by an undergraduate student, can be used to understand how is possible to tackle the Equation (1b). In fact, using the derivative of the inverse function, as in [18], or, alternatively, the theorem of implicit functions, as we are about to show in the Section 4.1, the student can learn how to build the power series representation for the Lambert function, taking advantage of Computer Algebra in a proactive way and not simply using it like a super pocket calculator. Our contribution suggests the development of a theme that leads the Advanced Calculus student to an organic reflection on multiple aspects of the subject, facilitating the solution of a certain class of problems connected with exponential-polynomial or logarithmic-polynomial equations. Finally, for completeness, it must be stated that the methods used for the numerical calculation of W are still the ongoing object of contemporary research [10]. We are strongly convinced that this approach helps form the right mindset in potential future applied scientists, who will use mathematical models in their profession.

Definition of W and its relationships with Calculus
To provide an idea on the nature of equations that can be solved using Lambert functions, we start with a working example consisting of a simple one dimensional Calculus problem. Example 1. Assuming that we are looking for the absolute minimum of f : R → R, f (x) = x 2 + e x . Since f is a coercive differentiable convex function, we infer "a priori" the existence of a unique minimiser, say x m which, of course, is such that 2x m + e x m = 0. At this point, in a traditional lesson, to identify this value the instructor refers to numerical methods, first and foremost the Newton-Raphson method. On the other hand, enquiring Mathematica ® we obtain: Based on this approach, the student learns that there is an explicit solution to the problem which, perhaps, they are not aware of. Moreover, Mathematica ® warns that, in any case, an issue exists which has to be understood. This is due to the fact that the Lambert function is, indeed, multivalued. This may, in our opinion, be one possible motivation for explaining how an explicit solution of this problem can be found.
We will limit our scope to the real line, having an undergraduate audience in mind. Lambert W functions are the multivalued inverse of : R → [−e −1 , +∞), (y) = ye y . Thus there are, indeed, two real Lambert functions, which are found solving, with respect to y, the transcendental equation provided that x ∈ [−e −1 , +∞). Thus W(x) is one of such solutions, obtained by inverting (y) in one of the We formally introduce the two real Lambert functions: The two branches of the Lambert function, W 0 and W −1 are constructed as follows: is the inverse function of (y), when restricted in the [−1, +∞) interval where it is strictly increasing, is the inverse function of (y), when restricted in the (−∞, −1) interval where it is strictly decreasing.
From now on, when we do not need to distinguish between the two branches, we write W(x) using the index, instead, when a certain statement holds true for only one of the two branches. By definition, for any y ∈ [−e −1 , +∞), then if ye y = x holds true, then: Notice that taking the logarithm of the first identity in (5), if x > 0, we get The Mathematica ® call for these functions is ProductLog[k,x] or, also, LambertW[k,x], taking k = 0 or k = −1; notice that the default is k = 0. Figure 3 shows both branches W −1 and W 0 . Although the basic derivation and integration formulas for W are well known, we report them for completeness of presentation. Since their proofs are provided in full detail in [7], we do not repeat them here. Theorem 1. The derivation rule for W is: Corollary 2. From the derivation rule (6) we infer that, as expected from its definition, when y → +∞ function W behaves like the logarithm. Indeed, we have Proof. To evaluate the (7) limit we can invoke De l'Hospital's Rule. The quotient of the derivatives is .
It is worth noting that the speed of convergence of (7) is extremely slow.

Remark 1.
Once again using De l'Hospital's Rule, it is possible to show that: Theorem 3. Discarding the constant of integration, the indefinite integration rule is: Knowing the primitive of W makes it possible to obtain more identities. For instance, in [13], the following relations for indefinite integrals are reported (again we omit the integration constant): Wikipedia provides two nice definite integrals connected with the probability integral and the Gamma function: We wish to add more integral relations here (the second was useful in [19], applied to an airplane trajectory model).
Proof. The (9a) and (9b) integrals are calculated using the substitution x = ue u = W −1 (u) so that dx = (1 + u)e u du, which leads to the u-integrals which, going back to the original variable, proves (9a) and (9b).
For the (9c) integral, we once again need a change of variable, in this case given by e x = z which transforms the integral on the right hand side of (9c) into the integral (8a), so that the identity (9c) is identified, returning to the original variable.

Transcendental equations
In this section, after some introductory working examples, we present our main results on transcendental polynomial-exponential and logarithmic-polynomial equations, solvable in terms of W. We will formulate general solution procedures that make it possible to symbolically treat a large class of such equations. Understanding the steps needed to solve these transcendental equations, which are hidden behind a blind use of Computer Algebra, undoubtedly represents a fundamental formative moment for students' mathematical maturity. Before illustrating the general results, let's take a look at a number of useful preparatory examples.
The first example appears naive at first glance. However, it is indeed quite informative on the use of W.

Example 2.
Consider, for x > 0, y > 0 the exponential equation The function f (y) = y y can be extended to [0, +∞) since f (y) → 1 as y → 0 + . We then observe that f (y) → +∞ as y → +∞ and that f (y) is C ∞ on (0, +∞). Since f (y) = y y (1 + ln y), the point (e −1 , e −1/e ) is the absolute minimum for f and therefore the Equation (10)  take the logarithm of both sides y y = x ⇐⇒ y ln y = ln x.
Now we change variable z = ln y so that we obtain, proceeding formally: In (11) we wrote W but, referring to Figure 4, if x ≥ 1, we have to choose W 0 while, if e −1/e < x < 1, there are two solutions of (10), since W −1 comes into play with W 0 . Going back to the original y = e z variable, we have the solution of (10)

Remark 2.
Using a similar argument, it is easy to see that the solution to the equation y y 2 = x is: for e −1/(2e) < x < 1, The following Propositions 2 and 3, and their corollaries provide the solution to a wide class of polynomial-exponential and polynomial-logarithmic equations. For this very reason, they are very useful in order to assist teaching activities.
Proof. The first statement follows graphing f (y) = αy − e y . It is immediate to verify that there is exactly one stationary point y = ln α, which is the absolute maximum since: Summing up, we see that for any y ∈ R we have f (y) ≤ f (ln α) = α (ln α − 1) . The plotting of f (y) is depicted in Figure 5 below. Thus, we have proved that if condition (12) holds, then Equation (13)  roots. To complete the proof, we show how they can both be represented in terms of the Lambert functions. First, we write (13) as follows: dividing both sides by α we get Then, taking the exponential and then multiplying both sides by −α we obtain: Introducing the temporary variable u = −α −1 e x (15) becomes: (16) is an equation in the "canonical" form (4) and, since the right hand side of (16) is negative, looking again at Figure 2, we see that there are two distinct solutions, expressed by W −1 and W 0 , if and only if the right hand side satisfies condition (17) below: but (17) is fulfilled due to the assumption that (12), thus (15) has two solutions, which we represent with the following compact notation: Next, recalling (13), we see that u = α −1 x − y and inserting in (18) we arrive at our thesis (14). Eventually, the last two statements are immediate.
Changing one sign in (13) we obtain a similar result but, this time, with only one real solution.
Proposition 3. Fix α > 0. Then equation: has a unique real root, given by: Proof. Here we deal with the function g(y) = αy + e y , therefore g (y) > 0 for any y and, being lim y→±∞ g(y) = ±∞, we infer that equation g(y) = x has a unique solution for any x ∈ R. To express this solution in terms of the W function, we proceed as before: divide both sides of (19) by α > 0 and then exponentiating we get: To reconstruct the structure of the (4) equation, we divide by α (21) and substituting u = α −1 e y we arrive at: Since the right hand side of (22) is positive, the solution is expressed in terms of W 0 . The Equation (20) follows repeating the final step of Proposition 2.

Remark 3. Given β = 0, if the Equation (13) is slightly modified into
and, similarly, Equation (19) is changed into it is immediate to replicate the proofs of Propositions 2 and 3 using the straightforward change of variables y 1 = βx, which produces equations of the same form consisting of (13) or (19) Once (25) is solved, it is immediate to arrive at the solutions of (23) and (24).

Remark 4 (Connection to the SIR model
where a, b, N, s(0) are positive parameters related to the model. Therefore using Proposition 2, we get: On page 413 of [20], it is stated that there is only one solution to the Equation (26) such that s(∞) < b/a . This is, indeed, confirmed by the fact that this property is verified by choosing W 0 in (27). Finally, referring to the existence conditions for solution (12) stated in Proposition 2, we stress that in the case of Equation (26)  Henceforth (26) has real solutions if the following condition is met: Propositions 2 and 3 are the natural starting points to deal with a larger family of equations, as illustrated below.

Remark 5.
Assuming α, β, x ∈ R, with αβ = 0 and b > 0, = 1. Then the equation: can be solved using Propositions 2 or 3, writing b y = e y ln b and changing the variable substituting y 1 = y ln b, thus transforming (28) into: Depending on the concordance of the signs of the two left hand side coefficients of the new unknown y 1 in (29), either Proposition 2 or Proposition 3 can be applied to solve (29) and therefore (28).
Finally, we can deal with polynomial logarithmic equations.

Remark 6.
Assuming α, β, b, x ∈ R, with αβ = 0 and b > 0, = 1. Then the equation: can be solved by first changing the variable substituting log b y = u and thus transforming (30) into: Then, to solve (31), we refer to Remark 5.

Inversion trough Implicit Function Theorem
After formalising, through a process of selective inversion, the procedure for recognising the Lambert function as the solution of a particular exponential-polynomial equation, it is useful to express the power series expansion of W. If the invertible function f : I → R, such that f (0) = 0 and f (0) = a 1 = 0, can be represented in a power series, for simplicity around the origin, as: The problem to express the link between the coefficients of the f Taylor series and those of f −1 naturally arises. This topic is not part of the standard undergraduate curricula and there is no doubt that it is an advanced level subject. To prove this fact, it is sufficient to note that it is reported in specialist monographs such as [21] paragraph 107 pp. 184-188 and [22] paragraph 7.3 pp. 129-133. the question was initially addressed by Bürmann and subsequently studied in a (rare) joint collaboration between A.M. Legendre and J.L. Lagrange [12]. Here, we recall Lagrange's formulation dating back to 1770 and published in the second volume of his "opera omnia" [23]. If f is a function given as in (32), that is f (0) = 0 and this zero is simple, then the solution of equation f (y) = x is given by x n n! d n−1 dy n−1 lim y→0 y f (y) n . (33) As an alternative to the (33) formula, it is possible, again following the early work of Lagrange, to compute higher derivatives of the inverse function as explained in the [24] and [25] articles. If g(x) denotes the inverse function of f (y) differentiating g ( f (y)) = y, we get g ( f (y)) f (y) = 1 and thus the derivatives of g follow in cascade: Notice that the relevant practical computations, expressed by (34), are effective only when a couple of if values (y, x), such that y = f (x) ⇐⇒ x = g(y), is known. It is interesting to note that in [3], the series inversion of the functions 1 (y) = ye −y is mentioned, which of course is closely related to our (y) = ye y and b(y) = y/(1 + y) 5 , as applications of Lagrange Brüman's formula (33), without however providing details. This, certainly, because the focus of the article does not concern this aspect, but rather to illustrate the hypergeometric nature of the solutions of the fifth degree equation, which, as is well known, cannot be expressed by radicals. In [18] this procedure is used to obtain the McLaurin series associated to the Lambert function. Instead, in this case to compute the higher derivatives of the inverse of a given function, and in particular for the Lambert function and b(y) proposed in [3], in order to express the solution of a particular quintic equation as a power series, we rely on the Implicit Function Theorem. In this case, we also make use of the fact that, thanks to Computer Algebra, we can easily and quickly overcome computational difficulties, taking advantage, moreover, of the opportunity to show to the more advanced student the use of the software. In support of the didactic effectiveness of the approach proposed here, we highlight two facts. First, proceeding using the successive derivatives of the inverse function as in (34), it is explicitly mandatory to calculate all the derivatives of the f function up to the desired order. Using, instead, the (33) formula, the calculations remain quite demanding. To give the reader an idea, we write down the first three addenda: x lim y 4y 2 f (y) 2 − y f (y) (y f (y) + 6 f (y)) + 2 f (y) 2 f 5 (y) .
Here, we not only have to evaluate higher derivatives of f , but we also have to deal with limits in undetermined form. For completeness of presentation and to define the notations, below we state the Implicit Function Theorem, in its two dimensional formulation.

Theorem 4.
Let Ω ⊂ R 2 be an open set and let F : Ω → R be a C 1 function. Suppose that there exists (x 0 , y 0 ) ∈ Ω such that F(x 0 , y 0 ) = 0 and F y (x 0 , y 0 ) = 0. Then, there exist δ , ε > 0 such that, for any x ∈ (x 0 − δ, and it subsequently holds that, for any x ∈ (x 0 − δ, x 0 + δ) : The inversion method based on the Implicit Functions Theorem is, by its nature, applicable to many situations, obviously respecting the fact that the implicit function which leads to the inversion is C ∞ . We will use the (35) formula of Theorem 4 in order to compute the higher order derivatives of ϕ, which represents the solution y of the equation F(x, y) = 0 parametrised by x. It is evident that, using (35), it is possible to obtain the derivatives computed in (34). In any case, as already stated, we prefer to operate with the (35) formula in particular situations (i.e. to provide the power series for W 0 in the next Section 4.1 and to describe the solution of the quintic equation x(y + 1) 5 − y = 0 in Section 4.2).

The McLaurin series for W 0
If we choose F(x, y) = (y) − x = ye y − x and (x 0 , y 0 ) = (0, 0), the hypotheses of Theorem 4 are fulfilled. At the same time, the conditions for the Lagrange inversion in terms of the (33) equation are also fulfilled, with these being (0) = 0, (0) = 1. The partial derivatives are easily computed as F y (x, y) = (1 + y)e y , F x (x, y) = −1, so we find the formula for the derivative of W 0 in an alternative way. Notice that, for ease of notation, here we write y(x) instead of W 0 (x) and, when there is not ambiguity, we shorthand this as y = y(x). Henceforth Observe that (36) does not contradict (6), since F (x, y(x)) = 0 means ye y = x. From (36), since y(0) = 0, it follows that y (0) = 1 and, using the chain rule, observing also that from (36) it follows that y is differentiable, we arrive at: The key point here, in starting a recursion procedure, is that we can take advantage of (36) and plugging it into (37) and representing the second derivative in terms of y only: (38) is "useful" because the sole knowledge of W 0 (0) = 0 allows us to compute W (2) 0 (0) = −2, as well as because we can differentiate it to obtain the third derivative and so on.
At this point, the problem of evaluating the derivatives of each order at the origin of y = W 0 has, in principle, been solved. To do this, however, very involved calculations are needed, which can be tackled using Computer Algebra. Again from (38), it follows that y (2) is differentiable and so on, for all the higher order derivatives. The procedure can be iterated and we can see that Mathematica ® has helped us a great deal to obtain y (3) and y (4) y (3) = 2y 2 + 8y + 9 e 3y (1 + y) 5 , y (4) = − 6y 3 + 36y 2 + 79y + 64 e 4y (1 + y) 7 .
This agrees with the known fact that the McLaurin series of W 0 (x) is The radius of convergence of (39) is 1/e: using the ratio test calling a n , the sequence of coefficients of x n on the right-hand side of (39), it follows that a n a n+1 = n n + 1 n−1 , and taking the limit for n → +∞. Observe that since 1/3 < 1/e, we can express the solution of (1b) in terms of an infinite series, in a similar way to its trigonometric counterpart (2), namely: (−1) n−1 n n−1 3 n n! .
In practice, the iterative procedure described here has allowed us to obtain the power series coefficients of the inverse of the (y) = ye y function without invoking the Lagrange-Bürmann formula (33). The didactic advantage is evident, since based on this approach the topic can be dealt with considerably in advance. The problem to provide a series representation of W, with an initial point different form the origin is an advanced one (such problem is tackled in [8]).

Approximate solution of a quintic equation
We conclude our article by applying the method to a quintic equation, which is approximated by extending the method seen for the Lambert function. The case of fifth degree equations is fascinating, since it is not possible, in generel, to express the solutions using radicals and therefore the approximate solution found assumes considerable importance. The complex issue of insolubility for radicals associated with a fifth degree equation is covered extensively in the [26] monograph. The quality of the approximation provided will be verified with an appropriate use of Mathematica ® . The solution of the quintic equation presented here, and more generally of the class of equations of degree greater or equal to five are expressed in terms of special functions, and therefore well beyond the undergraduate curricula. Here, relying on the Implicit Function Theorem, we provide an alternative way to compute the power series representation of the real solution of the equation where, as usual y is the unknown and x is a non-zero real parameter. In [3] it is stated that the solution to (40) is given by the superb formula: The reason behind the binomial formula (41) that expresses the (40) solution is explained by the theory of the differential resolvent. This theory was developed in the second half of the 19th century, starting from the contributions of [27][28][29][30][31][32][33][34][35]. The theory makes it possible to associate an algebraic equation with a linear differential equation with variable coefficients. The (41) formula can be obtained by integrating this differential equation using power series of hypergeometric nature: which are advanced topics, not treated at an undergraduate level, we refer to the classic monograph [36]. We also wish to mention the contribution of [37] where, not surprisingly, the power series is obtained using the Lagrange Bürmann formula. Clearly, the theory of the differential resolvent cannot be considered an undergraduate argument. However, the verification, at least empirical, of the fact that the power series (41) solves quintic Equation (40) is possible with the implicit derivation technique, already used for the series representation of W 0 .
We see that the computed values y(0), y (0), y (2) (0), y (3) (0), y (4) (0) are compliant with the b 0 , b 1 , b 2 , b 3 , b 4 sequence of coefficients of x n in the right-hand side of (41). The iterative nature of the procedure is easily implemented in Mathematica ® . This allows us to verify the statement for larger and larger values of the index n. Even if it is not a formal proof, the student is brought to verify in a concrete way the validity of the proposition and, hopefully, intrigued and motivated to look deeper into the formal aspects not typically covered by the first year curriculum.