+ coeff = one / (sqrt (four*(fromInteger k) - two)) :: a
+
+
+-- | A matrix containing 'big_N' functions indexed by their
+-- element/number. Each row in the matrix represents a finite element
+-- (i.e. an interval in the mesh). Within row @i@, column @j@ contains
+-- the @j@th 'big_N' basis function.
+--
+-- Any given 'big_N' will probably wind up in this matrix multiple
+-- times; the basis functions don't change depending on the
+-- interval. Only the /number/ of basis functions does. Computing
+-- them this way allows us to easily construct a lookup \"table\" of
+-- the proper dimensions.
+--
+-- The second example below relies on the fact that @big_N 3@ and
+-- @big_N 6@ expand to Legendre polynomials (2,4) and (5,7)
+-- respectively and so should be orthogonal over [-1,1].
+--
+-- Examples:
+--
+-- >>> import Data.Vector.Fixed ( N5, N6 )
+-- >>> import Integration.Gaussian ( gaussian )
+-- >>> import Linear.Matrix ( Col5, fromList )
+-- >>> import Naturals ( N19 )
+--
+-- >>> let p = fromList [[3],[3],[5],[4],[5]] :: Col5 Int
+-- >>> let mesh = undefined :: Col5 (Double,Double)
+-- >>> let params = Params mesh p :: Params N5 N5 N19 Double
+-- >>> let big_ns = big_N_matrix :: Mat N5 N6 (Double -> Double)
+-- >>> let n1 = big_ns !!! (1,0)
+-- >>> let n4 = big_ns !!! (4,0)
+-- >>> n1 1.5 == n4 1.5
+-- True
+-- >>> let n1 = big_ns !!! (1,3)
+-- >>> let n2 = big_ns !!! (2,4)
+-- >>> gaussian (\x -> (n1 x) * (n2 x)) < 1e-12
+-- True
+--
+big_N_matrix :: (Arity m, Arity n, Algebraic.C a, RealField.C a)
+ => Mat m n (a -> a)
+big_N_matrix =
+ construct lambda
+ where
+ lambda _ j x = big_N (toInteger j) x
+
+
+
+
+-- | Derivatives of the 'big_N's, that is, orthogonal basis functions
+-- over [-1,1]. The test case below comes from Sage where the
+-- orthogonality of the polynomials' derivatives can easily be
+-- tested. The indices are shifted by one so that k=0 is the first
+-- basis function.
+--
+-- Examples:
+--
+-- >>> import qualified Algebra.Absolute as Absolute ( abs )
+--
+-- >>> let expected = 11.5757525403319
+-- >>> let actual = big_N' 3 1.5 :: Double
+-- >>> Absolute.abs (actual - expected) < 1e-10
+-- True
+--
+big_N' :: forall a. (Algebraic.C a, RealField.C a) => Integer -> a -> a
+big_N' k x
+ | k < 0 = error "requested a negative basis function"
+ | k == 0 = negate ( one / (fromInteger 2))
+ | k == 1 = one / (fromInteger 2)
+ | otherwise = coeff * ( legendre k x )
+ where
+ two = fromInteger 2
+ coeff = sqrt ((two*(fromInteger k) + one) / two) :: a
+
+
+-- | The matrix of (N_i' * N_j') functions used in the integrand of
+-- the stiffness/mass matrices.
+big_N's_matrix :: (Arity m, Arity n, Algebraic.C a, RealField.C a)
+ => Mat m n (a -> a)
+big_N's_matrix =
+ construct lambda
+ where
+ lambda i j x = (big_N' (toInteger i) x) * (big_N' (toInteger j) x)
+
+
+-- | Compute the global load vector F.
+--
+-- Examples:
+--
+-- >>> import Data.Vector.Fixed ( N3, N4 )
+-- >>> import Linear.Matrix ( Col4, frobenius_norm, fromList )
+-- >>> import Naturals ( N7 )
+--
+-- >>> let big_A = const (1::Double)
+-- >>> let c x = sin x
+-- >>> let f x = x*(sin x)
+-- >>> let bdy = Left (Dirichlet (0,1::Double))
+-- >>> let pde = PDE big_A c f bdy
+--
+-- >>> let i1 = (0.0,1/3)
+-- >>> let i2 = (1/3,2/3)
+-- >>> let i3 = (2/3,4/5)
+-- >>> let i4 = (4/5,1.0)
+-- >>> let mesh = fromList [[i1], [i2], [i3], [i4]] :: Col4 (Double,Double)
+-- >>> let pvec = fromList [[2],[3],[2],[1]] :: Col4 Int
+-- >>> let params = Params mesh pvec :: Params N4 N3 N7 Double
+--
+-- >>> let f1 = [0.0418]
+-- >>> let f2 = [0.0805]
+-- >>> let f3 = [0.1007]
+-- >>> let f4 = [-0.0045]
+-- >>> let f5 = [-0.0332]
+-- >>> let f6 = [-0.0054]
+-- >>> let f7 = [-0.0267]
+-- >>> let expected = fromList [f1,f2,f3,f4,f5,f6,f7] :: Col N7 Double
+-- >>> let actual = big_F pde params
+-- >>> frobenius_norm (actual - expected) < 1e-4
+-- True
+--
+big_F :: forall m n l a.
+ (Arity l, Arity m, Arity n,
+ Algebraic.C a, RealField.C a, ToRational.C a)
+ => PDE a
+ -> Params m n l a
+ -> Col l a
+big_F pde params =
+ ifoldl2 accum zero (big_N_matrix :: Mat m (S n) (a -> a))
+ where
+ accum :: Int -> Int -> Col l a -> (a -> a) -> Col l a
+ accum i j prev_F this_N =
+ prev_F + this_F
+ where
+ two = fromInteger 2
+ (x1,x2) = (mesh params) !!! (i,0)
+ q = affine_inv (x1,x2)
+ integrand x = ((f pde) (q x)) * (this_N x)
+
+ -- The pointer matrix numbers from 1 so subtract one here to
+ -- get the right index.
+ global_idx = ((pointer params) !!! (i,j)) - 1
+ this_F = construct lambda
+ lambda k _ = if k == global_idx
+ then (gaussian integrand)*(x2 - x1) / two
+ else fromInteger 0
+
+
+
+big_K_elem :: forall m n l a b.
+ (Arity l, Arity m, Arity n,
+ Algebraic.C a, RealField.C a, ToRational.C a)
+ => PDE a
+ -> Params m n l a
+ -> Int
+ -> Int
+ -> Mat l l a
+ -> b
+ -> Mat l l a
+big_K_elem pde params _ k cur_K _ =
+ ifoldl2 accum cur_K (big_N's_matrix :: Mat m (S n) (a -> a))
+ where
+ accum :: Int -> Int -> Mat l l a -> (a -> a) -> Mat l l a
+ accum i j prev_K these_N's =
+ prev_K + this_K
+ where
+ two = fromInteger 2
+ (x1,x2) = (mesh params) !!! (k,0)
+ q = affine_inv (x1,x2)
+ integrand x = ((big_A pde) (q x)) * (these_N's x)
+ -- The pointer matrix numbers from 1 so subtract one here to
+ -- get the right index.
+ global_row_idx = ((pointer params) !!! (k,i)) - 1
+ global_col_idx = ((pointer params) !!! (k,j)) - 1
+ this_K = construct lambda
+ lambda v w = if v == global_row_idx && w == global_col_idx
+ then (two/(x2 - x1))* (gaussian integrand)
+ else fromInteger 0
+
+
+
+-- | Compute the \"big K\" stiffness matrix. There are three
+-- parameters needed for K, namely i,j,k so a fold over a matrix will
+-- not do. This little gimmick simulates a three-index fold by doing a
+-- two-index fold over a row of the proper dimensions.
+--
+-- Examples:
+--
+-- >>> import Data.Vector.Fixed ( N3, N4 )
+-- >>> import Linear.Matrix ( Col4, frobenius_norm, fromList )
+-- >>> import Naturals ( N7 )
+--
+-- >>> let big_A = const (1::Double)
+-- >>> let c x = sin x
+-- >>> let f x = x*(sin x)
+-- >>> let bdy = Left (Dirichlet (0,1::Double))
+-- >>> let pde = PDE big_A c f bdy
+--
+-- >>> let i1 = (0.0,1/3)
+-- >>> let i2 = (1/3,2/3)
+-- >>> let i3 = (2/3,4/5)
+-- >>> let i4 = (4/5,1.0)
+-- >>> let mesh = fromList [[i1], [i2], [i3], [i4]] :: Col4 (Double,Double)
+-- >>> let pvec = fromList [[2],[3],[2],[1]] :: Col4 Int
+-- >>> let params = Params mesh pvec :: Params N4 N3 N7 Double
+--
+-- >>> let k1 = [6, -3, 0, 0, 0, 0, 0] :: [Double]
+-- >>> let k2 = [-3, 10.5, -7.5, 0, 0, 0, 0] :: [Double]
+-- >>> let k3 = [0, -7.5, 12.5, 0, 0, 0, 0] :: [Double]
+-- >>> let k4 = [0, 0, 0, 6, 0, 0, 0] :: [Double]
+-- >>> let k5 = [0, 0, 0, 0, 6, 0, 0] :: [Double]
+-- >>> let k6 = [0, 0, 0, 0, 0, 6, 0] :: [Double]
+-- >>> let k7 = [0, 0, 0, 0, 0, 0, 15] :: [Double]
+-- >>> let expected = fromList [k1,k2,k3,k4,k5,k6,k7] :: Mat N7 N7 Double
+-- >>> let actual = big_K pde params
+-- >>> frobenius_norm (actual - expected) < 1e-10
+-- True
+--
+big_K :: forall m n l a.
+ (Arity l, Arity m, Arity n,
+ Algebraic.C a, RealField.C a, ToRational.C a)
+ => PDE a
+ -> Params m n l a
+ -> Mat l l a
+big_K pde params =
+ ifoldl2 (big_K_elem pde params) zero col_idxs
+ where
+ col_idxs = fromList [map fromInteger [0..]] :: Row m a
+