(33S) Legendre Polynomials

[Albert Chan]

(05-18-2022 02:59 AM)Eddie W. Shore Wrote:  Legendre Polynomials

The value of Legendre Polynomials can be calculated using a closed formula from Rodrigues' formula:

P_n(x) = Σ( comb(n, k) * comb(n+k, k) * ((x - 1)/2)^k, k= 0, n)

I think P_n as polynomial of (x-1)/2, is meant for symbolic calculations.

BTW, mpmath use this form: hyp2f1(-n, n+1, 1, (1-x)/2).
For integer n, both forms are equivalent (hyp2f1 version can work with arbitrary n)

Code:

function hyp2f1(a,b,c,z)
    local k, s, t = 1, 1, 1
    repeat
        t = t*a*b*z/(c*k)
        s = s + t
        a, b, c, k = a+1, b+1, c+1, k+1
    until s == s-t
    return s
end

function P_bad(n,x) return hyp2f1(-n,n+1,1,(1-x)/2) end

Numerical calculations with big n, coefs blow-up, with alternate signs, we have massive cancellations.
Below example, terms grow as big as ±1E48, before tappering down, making sum garbage.

P₁₀₀(0.25) = 1 - 3787.5 + 3585578.90625 - 1508034728.3203125 + ...

lua> P_bad(100, 0.25)
9.184297382421194e+031

Code:

function P(n,x) -- legendre recurrence relation formula
    local p0, p1 = 1, x
    for k=2,n do k0=k-1; k1=k+k0; p0, p1 = p1, (k1*p1*x-k0*p0)/k end
    return p1
end

Reccurence relation formula is faster, and *much* more accurate.



It may seems we have the same cancellation issues, but half the subtractions are really sums.
The other half, most of subtractions does not suffer catastrophic cancellations.

lua> P(100, 0.25) -- error = -2 ULP
0.07812465474298284