Third Order Newtons Method for Zernike Polynomial Zeros

where f (x), f ′ (x) and f ′′ (x) are the function and its first and second derivatives at the current best approximation x i [13,16,19]. For some classes of orthogonal polynomials, f ′′ /f ′ can be derived from f /f ′ [17,33,39], which means the update can be done to third order at essentially no additional numerical expense. If we divide the differential equation of the classical orthogonal polynomials, for example as tabulated in [1, 22.6][23],

(2)

Structure relations [24] relate the ratio f /f ′ to ratios at shifted indices n as tabulated for example in [1, 22.8],

The benefit is that the three-term recurrence equations, in the notation of [1, 22.7] (6)

. This is recursively inserted into the denominator of ( 5) to lower the index n until f 0 /f 1 is reached, which avoids problems with cancellation of digits. This work here implements this strategy for the Zernike polynomials, f = R m n , namely (i) fast calculation of f ′′ /f ′ from f /f ′ , (ii) calculation of f /f ′ from terminating continued fractions, both without evaluation of f or its derivatives via direct methods like Horner schemes.

2.1. Definition. We define Zernike radial polynomials in the unit Ball of dimension D in Noll’s nomenclature [29,31,20,37,5,35] for n ≥ 0, n-m = 0 (mod 2), m ≤ n as

x n-2s . (8) Remark 1. The normalization for D > 2 might be chosen differently [27].

Following the original notation, we will not put the upper (azimuth) index m in R m n -which is not a power-into parentheses. The normalization integral is (9)

Remark 2. The inversion of (8) assembles powers x i by sums of R m n (x), (i ≥ m, i -m even), (10) x i ≡ i n=m (mod 2)

By projection on the R m n basis these moments are

This terminating Saalschützian Hypergeometric Function has a closed form in terms of Γ-functions [22][34, (2.3.1.

3)]:

The round parenthesis with nonnegative lower indices in this equation are the finite products known as shifted factorials:

Much of this work is based on the representation as a terminating Gaussian Hypergeometric Function, the product of x m by a polynomial of degree n -m:

with the three “numerator” and “denominator” parameters

and argument

The initial terms of the power series (14) are

n are also Jacobi Polynomials [1, 15.4.6,22.5.42,22.5.1][3, 36]:

Remark 3. Gaussian integration rules for integrals

n changes sign over the integration interval. (i) (14) suggests to split R m n by assigning the factor x m to the weight such that a Gauss-Legendre integration for moments x D+m-1 is engaged and the wiggly remainder of R m n multiplied by f (x) is sampled over the abscissae. (ii) The representation (8) supports a hybrid term-wise method that adds the results of 1 + (n -m)/2 Gaussian integrations for moments x D-1+n-2s . The disadvantage here is the need for dense samples of f (x). (iii) There might be a workaround by developing rules for the lifted integrals

Elimination of F (a + ) establishes the format matching n-increments and n-decrements (14):

Substitution of ( 14) After insertion of these three formulas into ( 22)-( 24), the derivatives of

where ∼ = means the binomial factor and the argument list (a, b; c; z) of the hypergeometric function have not been written down explicitly.

Remark 4. These equations write the vector of derivatives of R m n as a lower triangular matrix with monomials of x multiplied by the vector of derivatives of F . Matrix inversion gives the reciprocal relations:

x -m-2 R ′′ ;

Since R m n (x) is a polynomial of order n, the (n + 1)st derivatives equal zero. Backward elimination of F and its derivatives with the aid of [1, 15.5.1] ( 28)

leads to the analog of ( 2),

This is one special case of differential equations that generate orthogonal functions [25], and could also be obtained by applying the derivatives of [1, 22.6.1] [10, 7] to (18). The derivative of this reaches out to the third derivatives, in which R ′′ is reduced to R and R ′ substituting the previous equation,

2.4. Zeros.

.

The analog of ( 5) is implemented by substituting [1, 15.

in the denominator. In lieu of ( 7) we find the continued fractions [12] (33)

which terminate in our cases since a is a negative integer and c = a + b. This already suffices to implement the standard Newton iteration, i.e., to approximate (1) by ∆x = -f (x)/f ′ (x). Division of (29

This is f ′′ (x)/f ′ (x) of the generic formula, and can be quickly computed from

of the lower order.

2.4.2. Initial Guesses. For n and m fixed, the strategy adopted here is to compute the (n -m)/2 distinct roots in (0, 1) starting with the smallest, then bootstrapping the others in naturally increasing order. An approximation to the smallest root is found by equating the first three terms in the square bracket of (

…(Full text truncated)…