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 where
  26
  27 import NumericPrelude
  28 import qualified Algebra.RealField as RealField
  29 import qualified Prelude as P
  30
  31
  32 -- | The @n@th Legendre polynomial in @x@ over [-1,1]. These are
  33 --   orthogonal over [-1,1] with the weight function w(x)=1.
  34 --
  35 --   Examples:
  36 --
  37 --   >>> let points = [0, 1, 2, 3, 4, 5] :: [Double]
  38 --
  39 --   The zeroth polynomial is the constant function p(x)=1.
  40 --
  41 --   >>> map (legendre 0) points
  42 --   [1.0,1.0,1.0,1.0,1.0,1.0]
  43 --
  44 --   The first polynomial is the identity:
  45 --
  46 --   >>> map (legendre 1) points
  47 --   [0.0,1.0,2.0,3.0,4.0,5.0]
  48 --
  49 --   The second polynomial should be (3x^2 - 1)/2:
  50 --
  51 --   >>> let actual = map (legendre 2) points
  52 --   >>> let f x = (3*x^2 - 1)/2
  53 --   >>> let expected = map f points
  54 --   >>> actual == expected
  55 --   True
  56 --
  57 --   The third polynomial should be (5x^3 - 3x)/2:
  58 --
  59 --   >>> let actual = map (legendre 3) points
  60 --   >>> let f x = (5*x^3 - 3*x)/2
  61 --   >>> let expected = map f points
  62 --   >>> actual == expected
  63 --   True
  64 --
  65 legendre :: (RealField.C a)
  66          => Integer -- ^ The degree (i.e. the number of) the polynomial you want.
  67          -> a -- ^ The dependent variable @x@.
  68          -> a
  69 legendre 0 _ = fromInteger 1
  70 legendre 1 x = x
  71 legendre n x =
  72   (c1*x * (legendre (n-1) x) - c2*(legendre (n-2) x)) / (fromInteger n)
  73   where
  74     c1 = fromInteger $ 2*n - 1
  75     c2 = fromInteger $ n - 1