Representing affine transform of Legendre polynomials

Question

I have a function defined as a set of weighted Legendre polynomials: $f(x)=\alpha_0 P_0(x) + \alpha_1 P_1(x) + \alpha_2 P_2(x) +\ldots$.

I have another function similarly defined with Legendre basis weights, but it's over a different range. I want to compare the two functions over the range where they overlap. This means I need to transform one function by stretching and shifting it to match the other sequence's range.

This is not too hard. I basically want to express my function with an affine transform applied to the input, and simplify the resulting polynomial down to basic Legendre basis weights. That is, finding the $\beta_i$ coefficients below.

$\begin{align*} f(Ax+B) &= \alpha_0 P_0(Ax+B) + \alpha_1 P_1(Ax+B) +\alpha_2 P_2(Ax+B) +\ldots\\ &= \beta_0 P_0(x) + \beta_1 P_1(x) + \beta_2 P_2(x) + \ldots \end{align*}$

This is quite well defined. I could in theory do it by expanding the sum into a simple giant polynomial, then collect terms and find the Legendre weights again.

So for example, for the second order term only,

$\begin{align*} P_2(Ax+B)&= {1\over2}(3(Ax+B)^2-1)\\ P_2(Ax+B)&= A^2(\frac12(3x^2-1)+ABx+{1\over2}(A^2+B^2)-\frac12\\ P_2(Ax+B)&= A^2 P_2(x) + ABP_1(x)+ {1\over 2}(A^2+B^2-1)P_0(x) \end{align*}$

This gets more and more ugly with higher orders.

My question: is there a nice recurrence relation or simple way to generate these transformed basis weights given A, B, and the set of weights $\alpha_i$, producing the resultant set of weights $\beta_i$? It feels like there should be a method that does not drop down to expanding the polynomial in full.

Answer 1

Well this is perhaps pretty much what you are doing, but it can be programmed quite simply. (Sorry about the lo-tech look, but I'm ASCII only).

The Legendre polynomials satisfy $P_{n+1} = d_n \cdot X\cdot P_n + f_n \cdot P_{n-1}$. If we define polynomials $p_n (x) = P_n(ax+b)$ then we have $p_n = \sum\limits_{j \leq n} \mathbf A_{n,j} P_j$ for a ( lower triangular) matrix $\mathbf A$, and the coefficient vector you seek will be $\mathbf A^T\dot \alpha$. Since the $p$'s also satisfy a recurrence relation, the rows of $\mathbf A$ can be computed recursively. I tried this and got:

A(0,0) = 1
A(1,0) = b    A(1,1) = a
n= 1..N-1
    A(n+1, 0) = p[n] * b * A( n , 0)
              - p[n] * a * (q[1]/p[1]) * A( n , 1)
              + q[n] * A( n-1, 0)

    j=1..n-1
        A(n+1, j) = p[n] * b * A( n, j) 
                  + p[n] * (a/p[j-1]) * A( n, j-1)
                  - p[n] * (a * q[j+1]/p[j+1]) * A( n, j+1)
                  + q[n] * A( n-1, j)

    A(n+1, n) = p[n] * b * A( n  , n)
              + p[n] * (a/p[n-1]) * A( n , n-1)

    A(n+1, n+1)= a*A( n , n)

which appears to work, though I cannot vouch for stability and suspect there has to be a better way.

Answer 2

The problem of converting from one orthogonal polynomial expansion to another has certainly been looked at before: Richard Hamming gave a conversion algorithm that turns a polynomial expressed in the monomial basis $x^i$ to a polynomial expressed in the Chebyshev basis $T_n(x)$, while Herbert Salzer generalized Hamming's conversion algorithm to convert from one orthogonal expansion to another. It is the Salzer algorithm which will be useful for your problem. (Actually, dmuir's proposal looks to be the Salzer algorithm, but slightly obscured.)

The theory behind Salzer's algorithm will take us too far afield to discuss (it involves the concept of so-called "mixed moments" between the old and new polynomial bases, and is tightly related to Chebyshev's algorithm which I mentioned here); suffice it to say that there exists a five term recurrence to yield the coefficients $\beta_j$ of $\sum\limits_{j=0}^n \beta_j P_j(x)$ from the coefficients $\alpha_j$ of $\sum\limits_{j=0}^n \alpha_j p_j^{(A,B)}(x)$. The formulae are too long for me to $\TeX$ out in full (see Salzer's paper for the formulae), but the implementation is not terribly complicated (and in fact you don't need a two-dimensional array, but only two one-dimensional arrays for temporary storage).

Here is a Mathematica implementation of your special case of Salzer's algorithm:

legcnv[cofs_?VectorQ, A_, B_] := 
  Module[{n = Length[cofs] - 1, c = B/A, dd, dp, d, i, j, k, v, w},
    dp = PadRight[{Last[cofs]}, n + 1];
    dd = PadRight[(2n - 1)Last[cofs]{B, A}/n + {cofs[[n]], 0}, n + 1];
    Do[
      i = n - k + 1;
      w = A(2i - 1)/i;
      v = dp[[1]]; d = dd[[1]]; dp[[1]] = d;
      dd[[1]] = cofs[[i]] + w(dd[[2]]/3 + c d) - i v/(i + 1);
      Do[
        v = dp[[j]]; d = dd[[j]]; dp[[j]] = d;
        dd[[j]] = 
          w(j dd[[j + 1]]/(2j + 1) + c d + (j - 1)dp[[j - 1]]/(2j - 3)) - 
            i v/(i + 1);
        , {j, 2, k - 1}];
      d = dd[[k]]; dp[[k]] = d;
      dd[[k]] = w((k - 1)dp[[k - 1]]/(2k - 3) + c d);
      dd[[k + 1]] = w k d/(2k - 1);
      , {k, 2, n}];
    dd]

(this should be easily translatable to any language).

Here's a symbolic test of Salzer's algorithm:

With[{n = 6}, Simplify[legcnv[u /@ Range[0, n], A, B].LegendreP[Range[0, n], x] == 
    (u /@ Range[0, n]).LegendreP[Range[0, n], A x + B]]]
True

and a numerical test:

legcnv[{4, 2, 1, 3, 1}, 3, -2]
{254, -576, 549, -297, 81}

that is, $4+2P_1(3x-2)+P_2(3x-2)+3P_3(3x-2)+P_4(3x-2)=254-576P_1(x)+549P_2(x)-297P_3(x)+81P_4(x)$