From: Michael Orlitzky <michael@orlitzky.com>
Date: Wed, 12 Feb 2014 05:00:27 +0000 (-0500)
Subject: Add big_K code to FEM.R1.
X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=2b6f975610bcc83e0c4a2fd26cd20fbcec06248d;p=numerical-analysis.git

Add big_K code to FEM.R1.
---

diff --git a/src/FEM/R1.hs b/src/FEM/R1.hs
index 02eab31..098dcf3 100644
--- a/src/FEM/R1.hs
+++ b/src/FEM/R1.hs
@@ -14,6 +14,17 @@
 --   if c(x) = 0, then it is assumed that the boundary conditions are
 --   Dirichlet.
 --
+--   The code creates a linear system,
+--
+--     (K + M)x = F
+--
+--   to be solved for @x@, where @K@ (\"big_K\") is the stiffness
+--   matrix, @M@ (\"big_M\") is the mass matrix, and @F@ (\"big_F\")
+--   is the load vector.
+--
+--   Warning: until PLU factorization is implemented, we can only
+--   solve the resulting system if it's positive definite!
+--
 module FEM.R1
 where
 
@@ -29,8 +40,10 @@ import Integration.Gaussian ( gaussian )
 import Linear.Matrix (
   Col,
   Mat(..),
+  Row,
   (!!!),
   construct,
+  fromList,
   ifoldl2,
   nrows )
 import Polynomials.Orthogonal ( legendre )
@@ -172,9 +185,10 @@ affine_inv (x1,x2) x =
     two = fromInteger 2
 
 
--- | Orthonormal basis functions over [-1,1]. The test case below
---   comes from Sage where the orthogonality of the polynomials'
---   derivatives can easily be tested.
+-- | Normalized integrals of orthogonal basis functions over
+--   n[-1,1]. The test case below comes from Sage where the
+--   orthogonality of the polynomials' derivatives can easily be
+--   tested.
 --
 --   Examples:
 --
@@ -215,7 +229,7 @@ big_N k x
 --
 --   Examples:
 --
---   >>> import Data.Vector.Fixed ( N5 )
+--   >>> import Data.Vector.Fixed ( N5, N6 )
 --   >>> import Integration.Gaussian ( gaussian )
 --   >>> import Linear.Matrix ( Col5, fromList )
 --   >>> import Naturals ( N19 )
@@ -223,7 +237,7 @@ big_N k x
 --   >>> 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 params
+--   >>> 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
@@ -233,17 +247,50 @@ big_N k x
 --   >>> gaussian (\x -> (n1 x) * (n2 x)) < 1e-12
 --   True
 --
-big_N_matrix :: (Arity m, Arity n, Arity l, Algebraic.C a, RealField.C a)
-             => Params m n l a
-             -> Mat m (S n) (a -> a)
-big_N_matrix params =
+big_N_matrix :: (Arity m, Arity n, Algebraic.C a, RealField.C a)
+             => Mat m n (a -> a)
+big_N_matrix =
   construct lambda
   where
-    maxdeg i = (max_degrees params) !!! (i,0)
-    lambda i j = if j > maxdeg i
-                 then (const $ fromInteger 0)
-                 else big_N (toInteger j)
+    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.
@@ -287,7 +334,7 @@ big_F :: forall m n l a.
       -> Params m n l a
       -> Col l a
 big_F pde params =
-  ifoldl2 accum zero (big_N_matrix 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 =
@@ -305,3 +352,86 @@ big_F pde params =
         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
+