1 {-# LANGUAGE ScopedTypeVariables #-}
3 -- | The polynomials over [a,b] form a basis for L_2[a,b]. But often
4 -- the \"obvious\" choice of a basis {1, x, x^2,...} is not
5 -- convenient, because its elements are not orthogonal.
7 -- There are several sets of polynomials that are known to be
8 -- orthogonal over various intervals. Here we provide formulae for a
9 -- few popular sets (bases).
11 -- First we present the simple formulae that generates polynomials
12 -- orthogonal over a specific interval; say, [-1, 1] (these are what
13 -- you'll find in a textbook). Then, we present a \"prime\" version
14 -- that produces polynomials orthogonal over an arbitrary interval
15 -- [a, b]. This is accomplished via an affine transform which shifts
16 -- and scales the variable but preserves orthogonality.
18 -- The choice of basis depends not only on the interval, but on the
19 -- weight function. The inner product used is not necessarily the
20 -- standard one in L_2, rather,
22 -- \<f, g\> = \int_{a}^{b} w(x)*f(x)*g(x) dx
24 -- where w(x) is some non-negative (or non-positive) weight function.
26 module Polynomials.Orthogonal (
31 import qualified Algebra.RealField as RealField ( C )
35 -- | The @n@th Legendre polynomial in @x@ over [-1,1]. These are
36 -- orthogonal over [-1,1] with the weight function w(x)=1.
40 -- >>> let points = [0, 1, 2, 3, 4, 5] :: [Double]
42 -- The zeroth polynomial is the constant function p(x)=1.
44 -- >>> map (legendre 0) points
45 -- [1.0,1.0,1.0,1.0,1.0,1.0]
47 -- The first polynomial is the identity:
49 -- >>> map (legendre 1) points
50 -- [0.0,1.0,2.0,3.0,4.0,5.0]
52 -- The second polynomial should be (3x^2 - 1)/2:
54 -- >>> let actual = map (legendre 2) points
55 -- >>> let f x = (3*x^2 - 1)/2
56 -- >>> let expected = map f points
57 -- >>> actual == expected
60 -- The third polynomial should be (5x^3 - 3x)/2:
62 -- >>> let actual = map (legendre 3) points
63 -- >>> let f x = (5*x^3 - 3*x)/2
64 -- >>> let expected = map f points
65 -- >>> actual == expected
68 legendre :: forall a. (RealField.C a)
69 => Integer -- ^ The degree (i.e. the number of) the polynomial you want.
70 -> a -- ^ The dependent variable @x@.
72 legendre 0 _ = fromInteger 1
75 (c1*x * (legendre (n-1) x) - c2*(legendre (n-2) x)) / (fromInteger n)
77 c1 = fromInteger $ 2*n - 1 :: a
78 c2 = fromInteger $ n - 1 :: a