Friday, July 19, 2013

Orthogonal Polynomials: Mathematica Script

A family of polynomials \[\{ \phi_{0}, \phi_1, ..., \phi_{n} \}\] (the subscript indicates the degree) are orthogonal with respect to a weighting function \(w(x)\), over an interval \([a,b]\), if \[\langle\phi_n, \phi_m\rangle = \int_{a}^{b} w(x) \phi_{n}(x) \phi_{m}(x) dx = \delta_{nm} c(n)\]
In other words, given a \(w(x)\) and an interval \([a, b]\) an orthogonal sequence of polynomials can be found via Gram-Schmidt orthogonalization of monomials.

This is done by a straightforward Mathematica script (link at the end). A possible use is:

P = OrthoPoly[Exp[-x], 0, 10, 3] // FullSimplify

This spits out the first three members of the family of orthogonal polynomials with \(w(x)=\exp(-x)\), over the interval [0,10].

One can use the zeros of these polynomials to obtain corresponding Gauss quadrature points. For example:

Solve[P[[3]] == 0, x] // N

yields {{x -> 0.180784}, {x -> 0.752396}}.

The code is available here.

No comments: