src/Polynomials/Orthogonal.hs

   1 {-# LANGUAGE ScopedTypeVariables #-}
   2
   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.
   6 --
   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).
  10 --
  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.
  17 --
  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,
  21 --
  22 --     \<f, g\> = \int_{a}^{b} w(x)*f(x)*g(x) dx
  23 --
  24 --   where w(x) is some non-negative (or non-positive) weight function.
  25 --
  26 module Polynomials.Orthogonal (
  27   legendre )
  28 where
  29
  30 import NumericPrelude
  31 import qualified Algebra.RealField as RealField ( C )
  32 import Prelude ()
  33
  34
  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.
  37 --
  38 --   Examples:
  39 --
  40 --   >>> let points = [0, 1, 2, 3, 4, 5] :: [Double]
  41 --
  42 --   The zeroth polynomial is the constant function p(x)=1.
  43 --
  44 --   >>> map (legendre 0) points
  45 --   [1.0,1.0,1.0,1.0,1.0,1.0]
  46 --
  47 --   The first polynomial is the identity:
  48 --
  49 --   >>> map (legendre 1) points
  50 --   [0.0,1.0,2.0,3.0,4.0,5.0]
  51 --
  52 --   The second polynomial should be (3x^2 - 1)/2:
  53 --
  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
  58 --   True
  59 --
  60 --   The third polynomial should be (5x^3 - 3x)/2:
  61 --
  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
  66 --   True
  67 --
  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@.
  71          -> a
  72 legendre 0 _ = fromInteger 1
  73 legendre 1 x = x
  74 legendre n x =
  75   (c1*x * (legendre (n-1) x) - c2*(legendre (n-2) x)) / (fromInteger n)
  76   where
  77     c1 = fromInteger $ 2*n - 1 :: a
  78     c2 = fromInteger $ n - 1 :: a