Thornahawk’s Numerical Analysis Notes

Most readers may be familiar with this particular definition of the natural logarithm:

The Lambert W function, named after J.H. Lambert who considered the solutions of a certain transcendental equation, is defined as follows:

Like the logarithm, Lambert’s function is multivalued, corresponding to the many possible complex solutions. Also like the logarithm, mathematicians have selected one of the function’s many possible branches to be the “principal branch”, sometimes denoted as . This particular branch is real for .

This function only very recently came into prominence compared to other, more popular special functions, due to the efforts of Corless, Knuth, and others.

A method for computing the function would involve solving for in the equation , and either of the iterative methods of Newton or Halley are suited for the task, due to the simple nature of the derivatives of the exponential function.

If Newton’s method is employed, the iteration becomes:

while the iteration using Halley’s method is:

The remaining problem is what to use as an initial approximation, or “seed” for either of the iterations to start with. Winitzki, in a paper on uniform approximations to transcendentals, presented the following expression:

where , and 4.6948, =0.8842, =0.9294, =0.5106, and =-1.213 are constants. This approximant matches the behavior of the Lambert W function at the points , , and at , and is accurate to around two significant figures.

Here, then, are TI 83+ programs for computing the principal branch of the Lambert W function. This one employs Newton’s method:

PROGRAM:LAMW1
√(2℮X+2)→V
2ln(1+.8842V)→Y
(Y–ln(1+.9294ln(1+.5106V))–1.213)/(1+(Y+.213ֿ¹)ֿ¹)→Y
If abs(Y+1)
Then
Repeat V=Y
Y→V
Y–(Y–X℮^(-Y))/(1+Y)→Y
ln(X)–ln(Y)→Y
End:End
Y

This one uses Halley’s method:

PROGRAM:LAMW2
√(2℮X+2)→V
2ln(1+.8842V)→Y
(Y–ln(1+.9294ln(1+.5106V))–1.213)/(1+(Y+.213ֿ¹)ֿ¹)→Y
If abs(Y+1)
Then
Repeat V=Y
Y→V
Y–X℮^(-Y)→W
Y–W/(1+Y–.5W(Y+2)/(1+Y))→Y
ln(X)–ln(Y)→Y
End:End
Y

Theoretically, the program using Halley’s method would be more efficient due to the faster convergence associated with the method; in practice, there is little difference in performance between the two for most complex arguments.

- thornahawk

References:

Corless, R.M., Gonnet, G.H., Hare, D.E.G., Jeffrey, D.J., and Knuth, D.E. “On the Lambert W Function.” Adv. Comput. Math. 5, 329-359, 1996.

Winitzki, S. “Uniform Approximations for Transcendental Functions.” In Computational Science and Its Applications – ICCSA 2003

The quadratic equation

is a very familiar sight to (most?) students and mathematicians alike. Finding the values of that solve the equation is a problem frequently encountered in applications.

One can solve the quadratic equation in either of two ways:

(the most frequently encountered form of the solution); or

(the form obtained by first dividing both sides of the equation by )

Most of the programs for solving the quadratic equation found at repositories like ticalc.org utilize more popularly the first form and to a certain extent the second form for their purpose.

Herein lies the problem of catastrophic cancellation. In brief, if , , or both are very small in magnitude compared to , one of the roots obtained through either form will be inaccurate since . The subtraction inherent in the computation of one of the roots will lead to a loss of significant figures.

Thus, one might want to somehow arrange things such that the operations involved will not have to depend on adding quantities of comparable magnitude and opposite sign.

The right way to utilize both forms of the quadratic equation is to use the following expression:

where one chooses the sign of the square root such that

and with that, the two solutions would be

If , , and are all expected to be real, can be expressed as

In TI-BASIC:

PROGRAM:QUADR
Prompt A,B,C
√(B²-4AC)→D
If B<0:-D→D
-.5(B+D)→Q
{Q/A,C/Q}→RTS


for the real version, and

PROGRAM:QUADC
Prompt A,B,C
√(B²-4AC)→D
If real(Dconj(B))<0
-D→D
-.5(B+D)→Q
{Q/A,C/Q}→RTS

for the complex version.

It is hoped that should one need to solve several quadratic equations, the pointers embodied in the two programs given be taken into account.

- thornahawk

References:

Press, W. H., et al. 1992, “Quadratic and Cubic Equations.” §5.6 in Numerical Recipes in FORTRAN: The Art of Scientific Computing, 2nd ed. (Cambridge, England: Cambridge University Press)

If you can see the quadratic formula and the power series for the exponential above this sentence, then the TeX rendering plugin has been successfully installed. \m/

The rendering plugin used for this blog was downloaded from here.

I will now be able to post what was meant to be placed in this blog. Keep watch.

- thornahawk

The name’s thornahawk. After three snags on manipulating Wordpress (mostly related to the problem of Permalink formatting), I hope that I’ll be able to smoothly add the meatier stuff from now on.

I conceived this blog as a repository of the notes I’ve gathered over the years of studying numerical methods, and the TI programs I’ve written to see how numerical algorithms work. I hope that what I’ll be writing may be both of interest and help to others.

I know this isn’t much of a welcome, but I try.

For now, contend with this post. I’ll still be studying the nuances of Wordpress before writing my thoughts down here.

Later.

- thornahawk