|
Note of
Numeric calculus |
Zeros of orthogonal polynomials
The NR method
A common method
for finding the zeros of orthogonal polynomials follows three basic steps
|
|
This task is performed with several approximate formulas: polynomials,
rational, or power function are usually adopted. |
|
This task is performed with recursive functions both for the polynomial
and the 1st derivative values. They do not use the polynomial coefficients
and are very fast end efficient |
|
|
This task is performed with the iterative use of the Newton-Raphson
formula. If the initial root approximation is sufficiently close the
convergence is very fast. But on the other hand it could converge to another
root if the initial estimation is not sufficiently accurate |
The key of
success is the first step. If the initial estimation of the root xi* is not sufficiently close
the algorithm will converge to another unwanted root. This may raise problems
when the degree of the polynomial increase because all the roots are more
closes each others.
The
following picture shows the zeros of Laguerre polynomials from 2 to 10 degree
As we can
see the distance between the roots increase at the right but decrease at left.
The distance between x1 and x2 becomes very small for
higher degree. This requires a good accuracy for the initial starting point of
the NR algorithm; otherwise it cannot find all the polynomial roots.
This
phenomenon explains why the estimation formulas are so complicated. For example
a good estimation formula for the generalized Laguerre polynomials known in the
literature ("Numerical recipes in fotran77 ...", Cambridge U Press,
1999) are the following
And its VB
implementation
If i = 1 Then
z(1) = (1 + m) * (3 + 0.92 * m) /
(1 + 2.4 * n + 1.8 * m)
ElseIf i = 2 Then
z(2) = z(1) + (15 + 6.25 * m) / (1 + 0.9 * m + 2.5 *
n)
Else
z(i) = z(i - 1) + ((1 + 2.55 * (i - 2)) / (1.9 * (i - 2)) + 1.26 * (i -
2) *
m / (1 + 3.5 * (i -
2)))* (z(i - 1) - z(i - 2)) / (1 + 0.3 * m)
Quite frightening,
isn't? To compensate this intrinsically complication, this formula is very fast
and accurate giving a global efficient method to calculate the zeros of the
generalized Laguerre polynomials
Fortunately
similar formulas exist also for other orthogonal polynomials. Following this
method
The DK method
An
alternative approche fo finding all zero of orthogonal polynomials use the so
called "
where z1,
z2, ... zr are a
subset of the polynomial roots. The iterative formula will converges to a polynomial
root surely different from the r-roots of the given subset. The thing works as
the r roots would be "extracted" from the polynomial p(z). This
avoids the convergence to these roots and explains the "implicit
deflating" term.
This method
does not require any accurate estimation starting point and it is also as fast
as the
This method
is valid for any orhtogonal polynomial family, (or better, for any polynomial
having all real roots).
Let's see
how it works
Attack strategy for orthogonal polynomials
Assume to
have to solve the 4 degree Laguerre polynomial.
We know
that the roots are all single:
0 < z1 < z2 < z3
< z4 < R
We can see
that R can be roughly estimated, for n £ 20 with the simple function
R @ 3.5×( n − 1)
Or with the
more accurate formula:
R @ 3.14*(n − 1)+0.018*(n − 1)2
At the
beginning, the subset is empty (r = 0) the DK formula comes back to the
If we
choose a starting point greater (but not too far) then the max root x4,
for example:
x0 = R = 3.5×(4-1) = 10.5,
The
iterative formula will converge monotonically to the max root x4
x4 @ 9.39507091230113
Note that
the starting point is not critical at all. The algorithm will work fine also
for x0= 12, 15, etc.
We have to
point out that we should not exceed with the starting point because orthogonal
polynomials raise sharply
When the
max root is found the DK formula changes in:
Now the new
starting point can be chosen near to x4, and greater of x3,
that it's the max root of the deflated polynomial.
This
situation is shown in the following plot.
Note that
we do not perform truly the deflating operation. The DK formula will simply
skip the x4 root.
A new starting
point satisfying the starting conditions x0 > x3 is:
x0 = x4 + h
Where
"h" is a little value (example 5% of x4). This is necessary because the DK formula
cannot started from a root itself (1/0 division error)
Substituting
we get the new estimate starting point
x0 @ 9.4 + 0.47 = 9.87
With this
starting point the DK formula will converge to the root x3,
obtaining
x3 @ 4.53662029692113
Again,
insert this value into the DK formula we have
And another
starting point
x0 @ 4.5 + 0.225 = 4.725
That will converge
to the root x2, obtaining
x2 @ 1.74576110115835
Again,
insert the x2 value into the DK formula we have
With the
starting point
x0 @ 1.7.5 + 0.05 = 1.785
That will
finally converge to the last root x1, obtaining
x1 @ 0.322547689619392
So, there
are the roots of 4th degree Laguerre polynomial
|
x4 = |
9.39507091230113 |
|
x3 = |
4.53662029692113 |
|
x2 = |
1.74576110115835 |
|
x1 = |
0.32254768961939 |
Now is
clear how this algorithm works and because the accuracy of all starting points
is not as important as in the NR method. We have only to bind at right the max
root and the DK algorithm will find all roots in sequence from the max to the
min.
|
|
This task is performed with
several approximate simple formulas. Accuracy of R is not required |
|
This task is performed with
recursive functions both for the polynomial and the 1st derivative values.
They do not use the polynomial coefficients and are very fast end efficient |
|
|
This task is performed with
the iterative use of the DK formula. If the starting point is greater than
the max root the convergence to the max root is guarateed. |
Incidentally
we note that all roots are already sorted. This is useful to save extra time
spent in the sorting step
A practical
example of this algorithm is implemented in OrthoPoly.xla
an Excel addin containing also macros
and functions for evaluating the
following polynomial families
|
Function Poly_ChebychevT(x, n) |
Chebychev
of 1st class |
|
Function Poly_ChebychevU(x, n) |
Chebychev
of 2nd class |
|
Function Poly_Gegenbauer(l, x, n) |
Gegenbauer |
|
Function Poly_Hermite(x, n) |
Hermite |
|
Function Poly_Jacobi(a, b, x, n) |
Jacobi |
|
Function Poly_Laguerre(x, n, m) |
Laguerre
generalized |
|
Function Poly_Legendre(x, n) |
Legendre |