Stieltjes polynomials

[Vick]

TL;DR Summary

orthogonal polynomials, Legendre polynomials

I am looking for a recurrence relation and/or defining expression for the Stieltjes polynomials with regard to the Legendre polynomials. I found an article about it here: Legendre-Stieltjes but they do not offer a formula.

For example a recurrence relation for the Gegenbauer polynomials is here: Gegenbauer under the second bullet point.

An example of a defining expression for the Jacobi polynomials is here: jacobi using hypergeometric function.

There is another article that deals with Stieltjes polynomials with regard to the Legendre function of the second kind here: Legendre function-Stieltjes found at the first set of Mathematica code. However as I don't understand the code, I don't really know what is the generating formula there.

 

[fresh_42]

Here are the Legendre polynomials
https://en.m.wikipedia.org/wiki/Associated_Legendre_polynomials
and here the Stieltjes modification
https://en.wikipedia.org/wiki/Stieltjes_polynomials

 

[Vick]

Here are the Legendre polynomials
https://en.m.wikipedia.org/wiki/Associated_Legendre_polynomials
and here the Stieltjes modification
https://en.wikipedia.org/wiki/Stieltjes_polynomials

Why associated Legendre polynomials?

For example the recurrence relation for the Legendre polynomials are:
$$
P_{0}(x) = 1
$$
$$
P_{1}(x) = x
$$
$$
nP_{n}(x) = (2n-1)xP_{n-1}(x)-(n-1)P_{n-2}(x)
$$

How do I modify them to obtain the Legendre-Stieltjes polynomials?

 

[Stephen Tashi]

I am looking for a recurrence relation and/or defining expression for the Stieltjes polynomials with regard to the Legendre polynomials. I found an article about it here: Legendre-Stieltjes but they do not offer a formula.

From some internet browsing this morning:

Articles such as http://www.mi.sanu.ac.rs/~gvm/radovi/f18-02.pdf state:

Orthogonal polynomials on the real line always satisfy a three-term recurrence relation.

However, the link you gave states:

The Legendre-Stieltjes polynomials do not satisfy three-term recurrence relations or have a particulary simple representation.

The link you gave defines the $(n+1)$th Legendre-Stieljes polynomial $E_{n+1}$ by:
$\int_{-1}^{1} E_{n+1} P_n(x) dx = 0$ for $n = 0,1,2,..n$ where $P_n(x)$ is the $n$th Legendre polynomial. [ Edit: correction $\int_{-1}^{1} E_{n+1}P_k(x) x^k = 0$ for $k = 0,1,2,...n$] As I visualize finding $E_k$ directly from that definition it would entail solving a system of $k$ equations in $k$ unknowns for the coefficients of $E_k$.

The link you gave cites the article https://www.ams.org/journals/mcom/1968-22-104/S0025-5718-68-99866-9/S0025-5718-68-99866-9.pdf. The remarks after eq. 3 indicate that solving the system of equations runs into problems of numerical accuracy. I don't know if those problems still exist in modern implementations of high precision arithmetic. If they do, then we could try to figure out the complicated looking methods described in that article following eq. 3.

 

[Vick]

The link you gave cites the article https://www.ams.org/journals/mcom/1968-22-104/S0025-5718-68-99866-9/S0025-5718-68-99866-9.pdf. The remarks after eq. 3 indicate that solving the system of equations runs into problems of numerical accuracy. I don't know if those problems still exist in modern implementations of high precision arithmetic. If they do, then we could try to figure out the complicated looking methods described in that article following eq. 3.

Yes, the article indeed says that there is neither a recurrence relation nor a simple representation. Can you please formulate the problem in math symbols and equations to find the system of k equations? Or else do you understand what they say their program does in order to find the Legendre-Stieltjes polynomials at the end? Something about constructor...it looks easy by the sound of it!

Otherwise there is a ready made computation in my last link which instead of finding the Stieltjes polynomials through the Legendre polynomials, it finds it through the Legendre function of the second kind. However their code is in mathematica. And I don't quite understand what the power series is doing?

The code is here:

SetAttributes[StieltjesE, Listable];
StieltjesE[n_Integer, x_] := Module[{xl},
     Set @@ Hold[StieltjesE[n, xl_],
                 Normal[Quiet[Series[4^(n - 1) Beta[n, n]/
                                     LegendreQ[n - 1, 0, 3, xl],
                                     {xl, ∞, 0}],
                              Series::sbyc]]];
     StieltjesE[n, x]]

Is it doing the following:

$$
E_{n}(x)= \sum_{k=0}^n 4^{k-1} {\frac {Beta(n,n)} {LegendreQ(k-1,0,3,x)}}
$$

 

[Stephen Tashi]

you please formulate the problem in math symbols and equations to find the system of k equations?

One of the principles of Algebra is that there is a nobility in suffering. I follow that principle only to a limited extent. I'll try an example showing one of the equations for one of the polynomials.

Consider the case $n = 3, k = 2$

Let $E_4(x) = c_4 x^4 + c_3 x^3 + c_2 x^2 + c_1 x + c_0$
where the $c_i$ are the unknown coefficients.

I didn't quote the correct definition for $E_n$ in my previous post. The correct definition is that:

$\int_{-1}^{1} E_{n+1}(x)P_k(x) x^k dx = 0$ for $k = 0,1,...n$So for $n =3$ :
$E_4(x) = E_{3+1}(x)$ satisifes $\int_{-1}^{1} E_4(x) P_k(x) x^k dx = 0$ for $k = 1,2,3$

To Derive the equation for the case $k = 2$ :

$P_2(x) = \frac{1}{2} (3 x^2 - 1)$

Let $f(x) = E_4(x) P_2(x) x^ 2 = (c_4 x^4 + c_3 x^3 + c_2 x^2 + c_1 x + c_0)(1/2)(3x^2 -1 ) x^2$

$= (1/2) \ ( \ c_4 (3x^6 - x^4) + c_3(3x^5 - x^3) + c_2(3x^4 - x^2) + c_1(3x^3 - x) + c_0(3x^2 - 1)\ ) (x^2)$

$= (1/2) \ (\ c_4( 3x^8 - x^6) + c_3(3x^7 - x^5) + c_2(3x^6 - x^4) + c_1(3x^5 - x^3) + c_0(3x^4 - x^2) \ )$

The antiderivative of $f(x)$ is:

$F(x) = \int f(x) dx = (1/2)\ (\ c_4( 3 \frac{x^{9}}{9} - \frac{x^7}{7}) + c_3(3 \frac{x^8}{8} - \frac{x^6}{6}) + c_2(3 \frac{x^7}{7} - \frac{x^5}{5}) + c_1(3 \frac{x^6}{6} - \frac{x^4}{4}) + c_0( 3 \frac{x^5}{5} - \frac{x^3}{3})\ )$

$\int_{-1}^1 f(x) dx = F(1) - F(-1)$

$= (1/2) ( c_4 (3/9 - 1/7 ) + c_3(3/8 - 1/6 ) + c_2(3/7 - 1/5 ) + c_1 (3/6 - 1/4 ) + c_0(3/5 - 1/3)$
$\ \ - (1/2) ( c_4( -3/9 + 1/7) + c_3( 3/8 - 1/6) + c_2(-3/7 + 1/5) + c_1 (3/6 - 1/4) + c_0(- 3/5 + 1/3) \ )$

$= (1/2) ( c_4( 6/9 - 2/7) + c_3(0) + c_2( 6/7 - 2/5) + c_1(0) + c_0(6/5 - 2/3)\ )$

Setting the above expression to zero gives the equation for $k = 3$ which is:

(6/9 - 2/7) c_4 + (6/7 - 2/5) c_2 + (6/5 - 2/3) c_0 = 0 ##

You should check my work. I don't have a reference to check it against.

Or else do you understand what they say their program does in order to find the Legendre-Stieltjes polynomials at the end? Something about constructor...it looks easy by the sound of it!

No, I don't. The documentation you linked does not give the code for the constructor. It only shows the syntax for calling it.

Otherwise there is a ready made computation in my last link which instead of finding the Stieltjes polynomials through the Legendre polynomials, it finds it through the Legendre function of the second kind. However their code is in mathematica. And I don't quite understand what the power series is doing?

I don't use Mathematica, so I don't know what that code does.

 

[Vick]

(6/9 - 2/7) c_4 + (6/7 - 2/5) c_2 + (6/5 - 2/3) c_0 = 0

For n=3 and k =1:
$c4(2/7)+c2(2/5)+c0(2/3)=0$
for n=3 and k =2:
$c4(6/9 - 2/7)+c2(6/7 - 2/5)+c0(6/5 - 2/3) = 0$
and for n=3 and k=3:
$c4(10/11 - 6/9)+c2(10/9 - 6/7)+c0(10/7 - 6/5) = 0$

But how do you find the coefficients if you are setting the 3 equations to equal 0?

BTW, I think you should multiply the equations for k=2 and 3 by (1/2) !

 

[Stephen Tashi]

But how do you find the coefficients if you are setting the 3 equations to equal 0?

It's a system of 3 simultaneous equations in 3 unknowns. Have you studied solving simultaneous linear equations?

BTW, I think you should multiply the equations for k=2 and 3 by (1/2) !

Multiplying or dividing both sides of an equation by a constant won't change its solutions.

The link you gave https://www.boost.org/doc/libs/1_66.../math_toolkit/sf_poly/legendre_stieltjes.html says:

E4(x) = P4(x) - 20P2(x)/27 + 14P0(x)/891

So $c_3 = c_1 = 0$ since the Legendre polynomials involved only have even powers.

We didn't get this conclusion merely from using the definition
$\int_{-1}^{1} E_4 P_k(x) x^k = 0$ For $k = 0,1, 2, 3$ (The link says $k = 0, 1,2...$ Other articles on the web begin with $k = 1$)

So the condition $\int_{-1}^{1} E_4 P_k(x) x^k$ is a property of $E_4$ and not a definition of $E_4$.

To understand a definition of $E_n$ we need to understand the first paragraphs of the link
https://www.ams.org/journals/mcom/1968-22-104/S0025-5718-68-99866-9/S0025-5718-68-99866-9.pdf

This has to do with Gaussian quadrature and I don't understand those paragraphs yet.

 

[Vick]

It's a system of 3 simultaneous equations in 3 unknowns. Have you studied solving simultaneous linear equations?

Yes I can solve any n simultaneous equations by matrix operations but as I said if all 3 equations = 0, all the coefficients will also equal 0.

∫1−1f(x)dx=F(1)−F(−1)∫−11f(x)dx=F(1)−F(−1)\int_{-1}^1 f(x) dx = F(1) - F(-1)

=(1/2)(c4(3/9−1/7)+c3(3/8−1/6)+c2(3/7−1/5)+c1(3/6−1/4)+c0(3/5−1/3)=(1/2)(c4(3/9−1/7)+c3(3/8−1/6)+c2(3/7−1/5)+c1(3/6−1/4)+c0(3/5−1/3)= (1/2) ( c_4 (3/9 - 1/7 ) + c_3(3/8 - 1/6 ) + c_2(3/7 - 1/5 ) + c_1 (3/6 - 1/4 ) + c_0(3/5 - 1/3)
−(1/2)(c4(−3/9+1/7)+c3(3/8−1/6)+c2(−3/7+1/5)+c1(3/6−1/4)+c0(−3/5+1/3) ) −(1/2)(c4(−3/9+1/7)+c3(3/8−1/6)+c2(−3/7+1/5)+c1(3/6−1/4)+c0(−3/5+1/3) )\ \ - (1/2) ( c_4( -3/9 + 1/7) + c_3( 3/8 - 1/6) + c_2(-3/7 + 1/5) + c_1 (3/6 - 1/4) + c_0(- 3/5 + 1/3) \ )

=(1/2)(c4(6/9−2/7)+c3(0)+c2(6/7−2/5)+c1(0)+c0(6/5−2/3) )=(1/2)(c4(6/9−2/7)+c3(0)+c2(6/7−2/5)+c1(0)+c0(6/5−2/3) ) = (1/2) ( c_4( 6/9 - 2/7) + c_3(0) + c_2( 6/7 - 2/5) + c_1(0) + c_0(6/5 - 2/3)\ )

(1/2) here is missing. Sympy says that k=2 = $c4* 4/21$
which is (6/9 -2/7) *0.5

 

[Stephen Tashi]

Yes I can solve any n simultaneous equations by matrix operations but as I said if all 3 equations = 0, all the coefficients will also equal 0.

Are you saying that is true only for these particular 3 equations? - or are you saying it is true for any 3 equations whose right hand side is zero?

 

[Vick]

Are you saying that is true only for these particular 3 equations? - or are you saying it is true for any 3 equations whose right hand side is zero?

Ok..there is homogeneous type...but it is really complicating matters, if I have to set the constants at 0!

What I meant before, normally we would have trivial solutions where all coefficients are zero as well!

 

[Stephen Tashi]

I quoted the definition incorrectly even in my edit!

The definition says $\int_{-1}^1 E_{n+1} P_n(x) x^k = 0$ for $k = 0,1,2,...,n$. There is no $k$ subscript on $P_n(x)$. So for the case $n = 3, k = 1$ we use $P_3(x)$. We don't use $P_1(x)$. Which $P$ did you use?

I realized this while transcribing quotations from https://www.ams.org/journals/mcom/1968-22-104/S0025-5718-68-99866-9/S0025-5718-68-99866-9.pdf. as a way of forcing myself to read them carefully.

2. The Extension of Quadrature Formulae. The basic reasoning behind the extension of quadrature formulae is as follows. Let an n-point formula be augmented by the addition of p abscissae and let $G_{n+p}(x)$ be the polynomial whose roots are the $n + p$ abscissae of the new quadrature formula. A general polynomial of degree $n + 2p -1$ can be expressed as

$(1)\$ $F_{n+2p-1}(x) = Q_{n+p-1}(x) + G_{n+p}(x) + \sum_{k=0}^{p+1}C_k x^k$

where $Q_{n+p-1}(x)$ is a general polynomial of degree $n+p-1$. This transformation of $F_{n+2p-1}(x)$ is possible since the number of unknown coefficients on the left- and right-hand sides of (1) is equal. $\ Q_{n+p-1}(x)$ can always be exactly integrated by a $(n + p)$-point formula and if $G_{n+p}(x)$ is such that

$(2)\ \$ $\int_{-1}^1 G_{n+p}(x)x^k dx = 0, \ k = 0,1,..., p-1$

then all of (1) can be exactly integrated by an $(n + p)$-point quadrature formula.Thus it should, in principle, be possible to derive formulae having $n+p$ abscissae and of degree $n + 2p - 1$.

2.1. The Extension of the Gauss Formulae. Kronrod [1] has considered the case $p = n + 1$ for the n-point Gauss formula.

For this choice of $p$ the polynomial $K{n+1}(x)$ whose $n + 1$ roots are the additional abscissae must satisfy, corresponding to (2),

$(3)\ \$ $\int_{-1}^{1} K_{n+1}(x)P_n(x)x^k dx = 0$ for $k = 0, 1,..., n$

, where $P_n(x)$ is the Legendre polynomial. Kronrod determines $K_{n+1}(x)$ and hence its zeros by substituting its polynomial expansion into (3) and solving the resulting triangular system of equations to find the polynomial coefficients. It is at this point that the numerical difficulties arise.

1. A. S. Kronrod, Nodes and Weights of Quadrature Formulas, English transl. from Russian, Consultants Bureau, New York, 1965. MR 32 #597

 

[Vick]

I quoted the definition incorrectly even in my edit!

The definition says $\int_{-1}^1 E_{n+1} P_n(x) x^k = 0$ for $k = 0,1,2,...,n$. There is no $k$ subscript on $P_n(x)$. So for the case $n = 3, k = 1$ we use $P_3(x)$. We don't use $P_1(x)$. Which $P$ did you use?

I realized this while transcribing quotations from https://www.ams.org/journals/mcom/1968-22-104/S0025-5718-68-99866-9/S0025-5718-68-99866-9.pdf. as a way of forcing myself to read them carefully.

I used the Legendre polynomials as the paper suggested, with reference to your help on how to go from the integration from -1 to 1 to a system of equations with constants = 0 and I solved it...

for k in range(1,r):
    b = legpoly(x,2*k-1)
    c = legpoly(x,n)
    print
    print "k = ",k
    print "------"
    for i in range(1,r+1):
        s = a+ str(m-i)
        f = (b*c*legpoly(x,2*i-1-q))
        g = risch_integrate(f,x)
        ii = g.subs(x,1)
        jj = g.subs(x, -1)
        if i + k < r:
            print s,0
        else:
            print s,ii-jj

The results are (taking m=4; n =3) : ...

k =  1
------
a3 0
a2 6/35
a1 8/63

k =  2
------
a3 2/7
a2 8/105
a1 4/77

As their article suggested to take the highest coefficient as unity, therefore we multiply a2 and a1 from k=2 by $7/2$

And using Gauss-Jordan Elimination to find the coefficients:

a = Matrix([[6.0/35,8.0/63,0],\
            [4.0/15,2.0/11,0]])

b = gjmat2(a,2)
pp.pprint(Matrix(b))

and the results are:

Matrix([
[1.0,  0.740740740740741, 0],
[  0, -0.015712682379349, 0]])

Thank you for helping me!