derivatives of Chebyshev polynomials

[Tamas Papp]

Hi,

Can someone point me to some Julia code that calculates a matrix for the
derivatives of the Chebyshev polynomials T_j, at given values, ie

d^k T_i(x_j) / dx^k for i=1,..n, j for some and vector x.

The Chebyshev polynomials themselves are very easy to calculate using
the recurrence relation, but derivatives are not. Alternatively, maybe
this can be extracted from ApproxFun but I have not found a way.

Best,

Tamas

[Spencer Lyon]

I have a function that will compute the derivative matrix operator that transforms coefficients of a Chebyshev interpolant on interval [a, b] of degree n (so n+1 total Chebyshev basis functions T_0, T_1, … T_n) into the coefficients of a degree n-1 chebyshev interpolant on [a, b] that is the exact derivative of the first interpolant. Maybe it could help you:

# derivative matrix
function der_matrix(deg::Int, a::Real=-1, b::Real=1)
    N = deg
    D = zeros(Float64, N, N+1)
    for i=1:N, j=1:N+1
        if i == 1
            if iseven(i + j)
                continue
            end
            D[i, j] = 2*(j-1)/(b-a)
        else
            if j < i || iseven(i+j)
                continue
            end
            D[i, j] = 4*(j-1)/(b-a)
        end
    end

    D
end

[Spencer Lyon]

Oh I forgot to say how to use it. If c are your coefficients for a degree n interpolant on [a, b] then the coefficients of the derivative of that interpolant are der_matrix(n, a, b) * c

[Alan Edelman]

there could be so many ways to do this

that it might be worth asking about the "use case"

the derivatives, of course, satisfy the same three term recurrence

but with different initial conditions

so it might be worth knowing if 

0.  whether this is worth optimizing a great deal or just a little bit or you just want something

 to do the job

1. the degree of the original polynomial was relatively small

2. whether one derivative in particular  or a run of many derivatives were desired

3. whether one or many points were desired

4. if the arguments were well inside or on [-1,1] 

[Paulo Jabardo]

There is a recurrence relation for the derivatives but it involves the Chebyshev polynomials of second kind.

Chebyshev polynomials of first kind:
T_0(x) = 1
T_1(x) = x
T_(n+1)(x) = 2x.T_n(x) - T_(n-1)(x)


Chebyshev polynomials of second kind:
U_0(x) = 1
U_1(x) = 2x
U_(n+1)(x) = 2x.U_n(x) - U_(n-1)(x)


Derivative:

dT_n/dx = n.U_(n-1)(x)




Source
http://en.wikipedia.org/wiki/Chebyshev_polynomials

NIST Handbook of mathematical functions:
http://dlmf.nist.gov/18.9



Paulo


PS I have a package https://github.com/pjabardo/Jacobi.jl
I will try to insert Chebyshev polynomials (and derivatives - very usefull for spectral methods) today.

[Alan Edelman]

If all you wanted was one derivative, i'd consider just an arccos and sines  formula

i gathered it was  multiple derivatives desired

[Tamas Papp]

Thanks to everyone for the replies.

I found that the excellent "Chebyshev and Fourier Spectral Methods" by J
P Boyd discusses computation of derivatives in one of the appendices
(A), suggesting 3 practical methods: (i) use the variable transformation
in terms of the cosine, (ii) use Gegenbauer polynomials, (iii) use the
derivative matrix transformation (as suggested by Spencer Lyon below).

My use case is roughly the following:

1. calculate the Chebyshev polynomials and the 1st (and possibly 2nd)
derivatives at a few (10-20) Chebyshev nodes,

2. use this to solve the functional equation (minimizing a residual,
Newton's method, etc),

3. evaluate the residual on a denser grid.

Since I only need to calculate the polynomials and their derivatives
once for 1., and once for 3., the cost is trivial.

Boyd gives a neat algorithm for method (i) above (Table 16.1, about 10
LOC, up to 4th derivative), so I went with that. IMO this is not worth
optimizing further. AFAIK method (iii) with the derivative matrix can be
slightly unstable numerically, there are quite a few papers on this.

Also, I decided to wrap up the code I wrote in a library, which I will
publish once I have experimented a bit with the interface and found it
satisfactory. I also coded up the standard routines for infinite and
semi-infinite domains to make solving PDEs less error prone.

Best,

Tamas