Content deleted Content added
TakuyaMurata (talk | contribs) m TakuyaMurata moved page Draft:Lentz's algorithm to Lentz's algorithm: likely notable |
m →Algorithm: fix "big K" notation |
||
| (37 intermediate revisions by 11 users not shown) | |||
Line 1:
In mathematics, '''Lentz's algorithm''' is an [[algorithm]] to evaluate [[generalized continued fraction|continued fraction]]s, and was originally devised to compute tables of spherical [[Bessel function]]s.<ref name=":0">{{Cite report|last=Lentz|first=W. J.|date=September 1973|title=A Method of Computing Spherical Bessel Functions of Complex Argument with Tables|url=https://apps.dtic.mil/sti/pdfs/AD0767223.pdf|publisher=Atmospheric Sciences Laboratory, US Army Electronics Command|___location=White Sands Missile Range, New Mexico|type=Research and Development Technical Report ECOM-5509 }}</ref><ref name="numerical-recipes-c++">{{cite book|title=Numerical Recipes in C++| pages=177–179|isbn= 0 521 75033 4}}</ref>
The version usually employed now is due to Thompson and Barnett.<ref name="Thompson-and-Barnett" />
== History ==
The idea was introduced in 1973 by William J. Lentz<ref name=":0" /> and was simplified by him in 1982.<ref>{{Cite book|last=J.|first=Lentz, W.|url=http://worldcat.org/oclc/227549426|title=A Simplification of Lentz's Algorithm.|date=August 1982|publisher=Defense Technical Information Center|oclc=227549426}}</ref> Lentz suggested that calculating ratios of spherical Bessel functions of complex arguments can be difficult. He developed a new continued fraction technique for calculating
== Initial work ==
This theory was initially motivated by Lentz's
== Algorithm ==
Lentz's algorithm is based on the [[Wallis-Euler relations]]. If
:<math>{f}_{0} = {b}_{0}</math>
:<math>{f}_{1} = {b}_{0} + \frac{{a}_{1}}{{b}_{1}}</math>
:<math>{f}_{2} = {b}_{0} + \frac{{a}_{1}}{{b}_{1} + \frac{{a}_{2}}{{b}_{2}}}</math>
:<math>{f}_{3} = {b}_{0} + \frac{{a}_{1}}{{b}_{1} + \frac{{a}_{2}}{{b}_{2} + \frac{{a}_{3}}{{b}_{3}}}}</math>
etc., or using the [[Generalized continued fraction#Notation|big-K notation]], if
:<math>{f}_{n} = {b}_{0} + \underset{j = 1}\overset{n}\operatorname{K}\frac{{a}_{j}}{{b}_{j}}</math>
is the <math>n</math>th convergent to <math>f</math> then
:<math>{f}_{n} = \frac{{A}_{n}}{{B}_{n}}</math>
where <math>{A}_{n}</math> and <math>{B}_{n}</math> are given by the Wallis-Euler recurrence relations
:<math>
\begin{align}
{A}_{-1} & = 1 & {B}_{-1} & = 0\\
{A}_{0} & = {b}_{0} & {B}_{0} & = 1\\
{A}_{n} & = {b}_{n} {A}_{n-1} + {a}_{n} {A}_{n-2} & {B}_{n} & = {b}_{n} {B}_{n-1} + {a}_{n} {B}_{n-2}
\end{align}
</math>
Lentz's method defines
:<math>{C}_{n} = \frac{{A}_{n}}{{A}_{n - 1}}</math>
:<math>{D}_{n} = \frac{{B}_{n - 1}}{{B}_{n}}</math>
so that the <math>n</math>th convergent is <math>{f}_{n} = {C}_{n} {D}_{n} {f}_{n - 1}</math> with <math>{f}_{0} = \frac{{A}_{0}}{{B}_{0}} = {b}_{0}</math> and uses the recurrence relations
:<math>
\begin{align}
{C}_{0} & = \frac{{A}_{0}}{{A}_{-1}} = {b}_{0} & {D}_{0} & = \frac{{B}_{-1}}{{B}_{0}} = 0\\
{C}_{n} & = {b}_{n} + \frac{{a}_{n}}{{C}_{n-1}} & {D}_{n} & = \frac{1}{{b}_{n} + {a}_{n} {D}_{n-1}}
\end{align}
</math>
When the product <math>{C}_{n} {D}_{n}</math> approaches unity with increasing <math>n</math>, it is hoped that <math>{f}_{n}</math> has converged to <math>f</math>.<ref name="Numerical-Recipes">{{Cite book |last1=Press |first1=W.H. |title=Numerical Recipes: The Art of Scientific Computing |last2=Teukolsky |first2=S.A. |last3=Vetterling |first3=W.T. |last4=Flannery |first4=B. P. |publisher=Cambridge University Press |year=2007 |edition=3rd |pages=207–208}}</ref>
Lentz's algorithm has the advantage of side-stepping an inconvenience of the Wallis-Euler relations, namely that the numerators <math>A_n</math> and denominators <math>B_n</math> are prone to grow or diminish very rapidly with increasing <math>n</math>. In direct numerical application of the Wallis-Euler relations, this means that <math>A_{n-1}</math>, <math>A_{n-2}</math>, <math>B_{n-1}</math>, <math>B_{n-2}</math> must be periodically checked and rescaled to avoid floating-point overflow or underflow.<ref name="Numerical-Recipes" />
== Thompson and Barnett modification ==
In Lentz's original algorithm, it can happen that <math>{C}_{n} = 0</math>, resulting in division by zero at the next step. The problem can be remedied simply by setting <math>{C}_{n} = \varepsilon</math> for some sufficiently small <math>\varepsilon</math>. This gives <math>{C}_{n + 1} = {b}_{n + 1} + \frac{{a}_{n + 1}}{\varepsilon} = \frac{{a}_{n + 1}}{\varepsilon}</math> to within floating-point precision, and the product <math>{C}_{n} {C}_{n + 1} = {a}_{n + 1}</math> irrespective of the precise value of ε. Accordingly, the value of <math>{f}_{0} = {C}_{0} = {b}_{0}</math> is also set to <math>\varepsilon</math> in the case of <math>{b}_{0} = 0</math>.
Similarly, if the denominator in <math>{D}_{n} = \frac{1}{{b}_{n} + {a}_{n} {D}_{n - 1}}</math> is zero, then setting <math>{D}_{n} = \frac{1}{\varepsilon}</math> for small enough <math>\varepsilon</math> gives <math>{D}_{n} {D}_{n + 1} = \frac{1}{{a}_{n + 1}}</math> irrespective of the value of <math>\varepsilon</math>.<ref name="Thompson-and-Barnett" /><ref name="Numerical-Recipes" />
== Applications ==
Lentz's algorithm was used widely in the late twentieth century. It was suggested that it doesn't have any rigorous analysis of error propagation. However, a few empirical tests suggest that it's
==References==
{{reflist}}
[[
| |||