src/Polynomials/Orthogonal.hs

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