[HP41] Lambert function RPN; question

[Thomas Klemm]

(12-26-2022 02:07 PM)Didier Lachieze Wrote:  Here is the original discussion from where this program is originated : Help solving equation.

It took me a while to figure out that I had used Halley's method with:

\(
\begin{align}
f(x) &= x e^x - a \\
f'(x) &= (1 + x) e^x \\
f''(x) &= (2 + x) e^x \\
\end{align}
\)

This is the formula used in the program:

\(
\begin{align}
x_{{n+1}}
&=x_{n}-{\frac {f(x_{n})}{f'(x_{n})-{\frac {f(x_{n})}{f'(x_{n})}}{\frac {f''(x_{n})}{2}}}} \\
\\
&=x_{n}-{\frac{x e^x - a}{(1 + x) e^x-\frac{(x e^x - a)(2 + x) e^x}{2(1 + x) e^x}}} \\
\\
&=x_{n}-{\frac{x e^x - a}{(1 + x) e^x-\frac{(x e^x - a)(2 + x)}{2(1 + x)}}} \\
\\
&=x_{n}-{\frac{x e^x - a}{(1 + x) e^x-\frac{(x e^x - a)(1 + \frac{1}{1 + x})}{2}}} \\
\end{align}
\)

But we could also use:

\(
\begin{align}
x_{n+1}
&=x_{n}-{\frac {f(x_{n})}{f'(x_{n})}}\left[1-{\frac {f(x_{n})}{f'(x_{n})}}\cdot \frac {f''(x_{n})}{2f'(x_{n})}\right]^{-1} \\
\\
&=x_{n}-{\frac {x e^x - a}{(1 + x) e^x}}\left[1-{\frac {x e^x - a}{(1 + x) e^x}}\cdot \frac {(2 + x) e^x}{2(1 + x) e^x}\right]^{-1} \\
\\
&=x_{n}+{\frac {a - x e^x}{e^x + x e^x}}\left[1+{\frac {a - x e^x}{e^x + x e^x}}\cdot \frac {2 + x}{2 + 2x}\right]^{-1} \\
\end{align}
\)

We end up with a slightly shorter program, here for the HP-42S:

Code:

00 { 41-Byte Prgm }
01▸LBL "W"
02 STO 00
03 0
04 STO 01
05▸LBL 00
06 RCL 00
07 RCL 01
08 ENTER
09 E↑X
10 ×
11 STO- ST Y
12 LASTX
13 +
14 ÷
15 RCL 01
16 2
17 RCL 01
18 +
19 +
20 LASTX
21 ÷
22 ÷
23 1
24 +
25 ÷
26 STO+ 01
27 ABS
28 RND
29 X≠0?
30 GTO 00
31 RCL 01
32 END

But it should work for the HP-41C as well.

Example

4 XEQ "W"

1.20216787320

So why Halley and not Newton, you might ask.
With a little more effort, we get cubic convergence, which is nice.