From: Michael Orlitzky 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 +