HP Forums

Full Version: Lambert W (HP-12C+)
You're currently viewing a stripped down version of our content. View the full version with proper formatting.

I have found a way to compute the the Lambert W function without using the solver, but only in the range (-1/e, 1]. However, the method allows the second branch to be computed all through the valid range, as shown below:

```--------------------
Lambert W (2nd branch), -1/e =< x < 0
01 CHS
02 1/x
03 STO 0
04 ENTER
05 ENTER
06 ENTER
07 STO 1
08 LN
09 *
10 RCL 1
11 x<>y
12 -
13 x=0
14 GTO 19
15 CLx
16 RCL 0
17 LST x
18 GTO 07
19 LST x
20 RCL 0
21 /
22 CHS
23 GTO 00
12C+            actual
-0.367879441 R/S ->  -1.000000000     -1.000030530
-0.367       R/S ->  -1.070791879     -1.070791887
-0.3         R/S ->  -1.781337023     -1.781337023
-0.000000001 R/S -> -23.897019580    -23.897019585
(actual results obtained with Egan Ford's LAMBERTWB function in his special functions library - See message #6 in
W-1(-1/x) = -(...x*ln(x*ln(x*ln(x*ln(x))))...)/x     (yet to be proved)
--------------------
```

The principal branch, at least in the range (-1/e, 1], can be computed with the following program:

```
--------------------
Lambert W (principal branch), -1/e > x >= 1
24 STO 0
25 ENTER
26 ENTER
27 ENTER
28 STO 1
29 ex
30 /
31 RCL 1
32 x<>y
33 -
34 x=0
35 GTO 40
36 CLx
37 RCL 0
38 LST x
39 GTO 28
40 LSX x
41 GTO 00
12C+           actual
0.999999999 GTO 24 R/S ->  0.567143290     0.567143290
0.5         GTO 24 R/S ->  0.351733711     0.351733711
0.000000000 GTO 24 R/S ->  0.000000000     0.000000000
-0.3         GTO 24 R/S -> -0.489402227    -0.489402227
-0.367       GTO 24 R/S -> -0.932399183    -0.932399185
(occasionally the program will loop for certain arguments, as for 1, due to rounding errors - in this case stop the program, the answer should be in register 1)
--------------------
```

Except for arguments whose convergences are very slow, these should work on the 12C as well. Also, by filling the stack with the desired argument and iterating two keystrokes until convergence is reached, this should work even on a non-programmable RPN calculator. The method may not be so practical but is somewhat interesting.

Gerson.

Edited: 8 Nov 2009, 7:12 a.m.

Gerson,

Thanks for the program and for the reference. I have often come across the Lambert function, but did not stop to read about it. Your program and the Wikipedia reference gave me the info I need to implement the Lambert function (if and when I feel like doing that). The approximations are interesting.

Namir

Hello Namir,

Quote:
The approximations are interesting.

Actually these are supposed to be exact solutions (except for accumulated rounding errors, mainly at the extremes of the convergence intervals). I've just found this paper which sheds some light on the subject (Iteration, pages 18 and 19).

The programs have been rewritten:

```-------------
01 CHS               20 STO 0
02 1/x               21 STO 1
03 STO 0             22 4
04 STO 1             23 EEX
05 LN                24 CHS
06 RCL 0             25 1
07 *                 26 0
08 RCL 1             27 STO 2
09 x<>y              28 RCL 0
10 -                 29 RCL 1
11 x=0               30 e^x
12 GTO 15            31 /
13 LST x             32 RCL 1
14 GTO 04            33 x<>y
15 LST x             34 STO 1
16 RCL 0             35 -
17 /                 36 ENTER
18 CHS               37 *
19 GTO 00            38 SQRT
39 RCL 2
40 x<=y
41 GTO 28
42 RCL 1
43 GTO 00
actual
-0.3678794411         R/S ->  -1.0000000000  (-1.0000197080)
1.0000000000 GTO 20  R/S ->   0.5671432905  ( 0.5671432904)
-------------
```

Regards,

Gerson.

Edited: 8 Nov 2009, 1:54 p.m.

I was playing with that a little.

It seems like

y(n+1)=x/exp(y(n))

converges to W(x) up to about 2.7 (probably supposed to be e). If you use

y(n+1)=(x/exp(y(n))+y(n))/2

it converges much faster in the previous range (because the oscillations are suppressed)(though it doesn't necessarily speed the convergence for negative argument), and also extends the range substantially, to about 60 (probably supposed to be 3*exp(3)).

y(n+1)=(x/exp(y(n))+2*y(n))/3

converges all the way to 5*exp(5)~742.

It's probably possible to use this or a similar idea to get an arbitrarily wide range of convergence.

In fact

y(n+1) = (x/exp(y(n)) + (y(n))^2)/(1+y(n))

seems to converge over the entire range. You just need a good seed for y(0). Maybe x for x < 1, log(x) for x>1. Or just the natural log of x+1 for all x.

Edited: 10 Nov 2009, 12:26 p.m.

Quote:
In fact

y(n+1) = (x/exp(y(n)) + (y(n))^2)/(1+y(n))

seems to converge over the entire range. You just need a good seed for y(0). Maybe x for x < 1, log(x) for x>1. Or just the natural log of x+1 for all x.

That's great! A good estimate for x>1 might be ln-ln(ln(x)). Perhaps your first expression and another initial estimate should be used for x<1.01 so that only a few iterations are needed in that range.

Some examples (Excel 2007)

```	yn+1 = (x/exp(yn) + yn2)/(1+yn)
x = 	1.001
n	yn
0	6.90925457110691
1	6.03581510741710
2	5.17828530065102
3	4.34105591640598
4	3.53072548897921
5	2.75791050589333
6	2.04091021461044
7	1.41252273893777
8	0.92807090671205
9	0.65196072145989
10	0.57301422876145
11	0.56752983893071
12	0.56750507994374
13	0.56750507944168
14	0.56750507944168
x = 	2
0	1.05966010114161
1	0.88171184629193
2	0.85324579958077
3	0.85260581752487
4	0.85260550201380
5	0.85260550201373
6	0.85260550201373
x = 	123456789012345
0	28.96730703740500
1	29.08342041664470
2	29.07698748392480
3	29.07696600979440
4	29.07696600955620
5	29.07696600955620
x = 	2.E+100
0	225.50944807065200
1	225.53347447862600
2	225.53318919246000
3	225.53318915157900
4	225.53318915157900
------------------
Formulas
x = 	2E+100
n	yn
0	=LN(D5)-LN(LN(D5))
1	=(\$D\$5/EXP(D9)+D9^2)/(1+D9)
2	=(\$D\$5/EXP(D10)+D10^2)/(1+D10)
------------------
```

Regards,

Gerson.

Edited: 10 Nov 2009, 4:38 p.m.

I notice y(n+1)=log(x/y(n)) converges to the second branch quickly if x is very small, but converges very, very slowly if y~-1. I thought a bit about speeding up the convergence, but am not quite there.

Incidentally, look at the end of this paper.

Quote:
Finally, in an eloge near the beginning of the two volumes of Lambert's mathematical papers (published in Zurich in 1946 and 1948) we find that Lambert had dreams of building a machine to do symbolic calculation (while Pascal's machine merely did arithmetic). Thus his wishes anticipated those of Lady Ada Lovelace by roughly one century, and the actuality of symbolic computation systems by roughly two.

I'd think Lambert would have to be impressed by the HP50g.

Okay, I think I have a good improvement to

y(n+1)=log(x/y)

Try

y(n+1)=(log(x/y)+1)*y/(1+y)

And maybe y(0)=log(-x)

Small x test

x=-.000001

Converges to -16.6265089014 in only three iterations. Takes 10 iterations for the other method.

Medium x test

x=-0.18

Converges to -2.7127679118 in only four iterations. Takes about 27 iterations to converge with other method.

Large x test

x=-0.36787944

Converges to -1.00007980968 in 17 iterations. Using qbasic and 100,000 iterations with the other method, has only converged to -1.00007975...

Edited: 10 Nov 2009, 11:17 p.m.

D'oh! That's just Newton's method to solve y-log(x/y)=0. Oh, well.

Quote:
y(n+1) = (x/exp(y(n)) + (y(n))^2)/(1+y(n))

The HP-12C program below uses your formula and y(0)=ln(x+2)-ln(ln(x+2)) for all x. It's very fast on the 12C+. Running times are acceptable on the original 12C.

```Lambert W function (principal branch)
01 STO 0	16 1		31 +
02 2		17 0		32 x<>y
03 +		18 x<>y		33 /
04 LN		19 y^X		34 RCL 1
05 ENTER	20 STO 2	35 x<>y
06 LN		21 1		36 STO 1
07 -		22 RCL 1	37 -
08 STO 1	23 +		38 ENTER
09 LN		24 LST x	39 *
10 1		25 ENTER	40 SQRT
11 0		26 * 		41 RCL 2
12 LN		27 RCL 0	42 x<=y
13 /		28 LST x	43 GTO 21
14 9		29 e^x 		44 RCL 1
15 -		30 /		45 GTO 00
12C+           actual    time (12C)
9.999999e99 R/S ->  224.8431063     224.8431064      09.5 s
99.999999990 R/S ->    3.385630139     3.385630140    13.9 s
2.000000000 R/S ->    0.852605502     0.852605502    13.7 s
1.000000000 R/S ->    0.567143290     0.567143290    16.4 s
0.5         R/S ->    0.351733711     0.351733711    15.9 s
0.000000000 R/S ->    5.765480-25     0.000000000    17.2 s
-0.3         R/S ->   -0.489402227    -0.489402227    22.4 s
-0.367       R/S ->   -0.932399184    -0.932399185    30.9 s
-0.000100000 R/S ->   -0.000100010    -0.000100010    17.8 s
```

Edited to include running times on the original HP-12C

Edited: 12 Nov 2009, 7:14 a.m.

I think this forum is the first place I've ever heard of the Lambert W. I've been reading up on it a bit. There's some interesting stuff about it.

Quote:
I think this forum is the first place I've ever heard of the Lambert W.

That's where I first heard of it too. I think it was Valentin Albillo who first mentioned it here.

I like its being capable of providing analytical solutions to some problems that would otherwise require a numerical solver, like the ridget hurler example in the HP-15C manual:

h = 5000(1 - e - t/20) - 200t

The problem consists in solving the equation for t when h = 0. If Lambert W were used, the solution would be simply

t = 20W(-5/4 e-5/4) + 25

Or t = 9.284255090, when evaluated with the HP-12C program above.

Edit to correct a typo.

Edited: 21 Nov 2009, 9:25 a.m.