From 1e988aee439f5352c25a4de03d7d18a58ae1d403 Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Sat, 15 Feb 2014 21:27:34 -0500 Subject: [PATCH] Add 'coefficients' to FEM.R1 which solves the desired system. --- src/FEM/R1.hs | 231 +++++++++++++++++++++++++++++++++++++++----------- 1 file changed, 182 insertions(+), 49 deletions(-) diff --git a/src/FEM/R1.hs b/src/FEM/R1.hs index 098dcf3..922469f 100644 --- a/src/FEM/R1.hs +++ b/src/FEM/R1.hs @@ -45,7 +45,9 @@ import Linear.Matrix ( construct, fromList, ifoldl2, - nrows ) + nrows, + set_idx ) +import Linear.System ( solve_positive_definite ) import Polynomials.Orthogonal ( legendre ) -- | Dirichlet boundary conditions. Since u(a)=u(b)=0 are fixed, @@ -185,6 +187,8 @@ affine_inv (x1,x2) x = two = fromInteger 2 +-- * Load vector + -- | 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 @@ -255,51 +259,22 @@ big_N_matrix = 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) +-- | The matrix of (N_i * N_j) functions used in the integrand of +-- the mass matrices. +big_Ns_matrix :: (Arity m, Arity n, Algebraic.C a, RealField.C a) => Mat m n (a -> a) -big_N's_matrix = +big_Ns_matrix = construct lambda where - lambda i j x = (big_N' (toInteger i) x) * (big_N' (toInteger j) x) + 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 ) +-- >>> import Naturals ( N3, N4, N7 ) -- -- >>> let big_A = const (1::Double) -- >>> let c x = sin x @@ -347,12 +322,47 @@ big_F pde params = -- 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 + k = ((pointer params) !!! (i,j)) - 1 + integral = (gaussian integrand)*(x2 - x1) / two + this_F = set_idx zero (k,0) integral + + +-- * Stiffness matrix + +-- | 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 matrix. +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) big_K_elem :: forall m n l a b. @@ -378,12 +388,10 @@ big_K_elem pde params _ k cur_K _ = 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 + row_idx = ((pointer params) !!! (k,i)) - 1 + col_idx = ((pointer params) !!! (k,j)) - 1 + integral = (two/(x2 - x1))* (gaussian integrand) + this_K = set_idx zero (row_idx, col_idx) integral @@ -394,9 +402,8 @@ big_K_elem pde params _ k cur_K _ = -- -- Examples: -- --- >>> import Data.Vector.Fixed ( N3, N4 ) -- >>> import Linear.Matrix ( Col4, frobenius_norm, fromList ) --- >>> import Naturals ( N7 ) +-- >>> import Naturals ( N3, N4, N7 ) -- -- >>> let big_A = const (1::Double) -- >>> let c x = sin x @@ -435,3 +442,129 @@ big_K pde params = where col_idxs = fromList [map fromInteger [0..]] :: Row m a + +-- * Mass matrix + +big_M_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_M_elem pde params _ k cur_M _ = + ifoldl2 accum cur_M (big_Ns_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_M these_Ns = + prev_M + this_M + where + two = fromInteger 2 + (x1,x2) = (mesh params) !!! (k,0) + q = affine_inv (x1,x2) + integrand x = ((c pde) (q x)) * (these_Ns x) + -- The pointer matrix numbers from 1 so subtract one here to + -- get the right index. + row_idx = ((pointer params) !!! (k,i)) - 1 + col_idx = ((pointer params) !!! (k,j)) - 1 + integral = (x2 - x1)*(gaussian integrand) / two + this_M = set_idx zero (row_idx, col_idx) integral + + +-- | Compute the \"big M\" mass matrix. There are three +-- parameters needed for M, 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 Linear.Matrix ( Col4, frobenius_norm, fromList ) +-- >>> import Naturals ( N3, N4, 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 m1 = [0.0723,0.0266,0,-0.0135,-0.0305,0.0058,0] :: [Double] +-- >>> let m2 = [0.0266,0.0897,0.0149,0,-0.0345,-0.0109,-0.0179] :: [Double] +-- >>> let m3 = [0,0.0149,0.0809,0,0,0,-0.0185] :: [Double] +-- >>> let m4 = [-0.0135,0,0,0.0110,0,0,0] :: [Double] +-- >>> let m5 = [-0.0305,-0.0345,0,0,0.0319,0.0018,0] :: [Double] +-- >>> let m6 = [0.0058,-0.0109,0,0,0.0018,0.0076,0] :: [Double] +-- >>> let m7 = [0,-0.0179,-0.0185,0,0,0,0.0178] :: [Double] +-- +-- >>> let expected = fromList [m1,m2,m3,m4,m5,m6,m7] :: Mat N7 N7 Double +-- >>> let actual = big_M pde params +-- >>> frobenius_norm (actual - expected) < 1e-3 +-- True +-- +big_M :: 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_M pde params = + ifoldl2 (big_M_elem pde params) zero col_idxs + where + col_idxs = fromList [map fromInteger [0..]] :: Row m a + + + +-- | Determine the coefficient vector @x@ from the system @(K + M)x = F@. +-- +-- Examples: +-- +-- >>> import Linear.Matrix ( Col4, Col7, frobenius_norm, fromList ) +-- >>> import Naturals ( N3, N4, 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 c1 = [0.02366220347687] :: [Double] +-- >>> let c2 = [0.03431630082636] :: [Double] +-- >>> let c3 = [0.02841800893264] :: [Double] +-- >>> let c4 = [-0.00069489654996] :: [Double] +-- >>> let c5 = [-0.00518637005151] :: [Double] +-- >>> let c6 = [-0.00085028505337] :: [Double] +-- >>> let c7 = [-0.00170478210110] :: [Double] +-- >>> let expected = fromList [c1,c2,c3,c4,c5,c6,c7] :: Col7 Double +-- >>> let actual = coefficients pde params +-- >>> frobenius_norm (actual - expected) < 1e-8 +-- True +-- +coefficients :: forall m n l a. + (Arity m, Arity n, Arity l, + Algebraic.C a, Eq a, RealField.C a, ToRational.C a) + => PDE a + -> Params m n (S l) a + -> Col (S l) a +coefficients pde params = + solve_positive_definite matrix b + where + matrix = (big_K pde params) + (big_M pde params) + b = big_F pde params -- 2.43.2