Row,
(!!!),
construct,
+ element_sum2,
fromList,
ifoldl2,
- nrows )
+ map2,
+ nrows,
+ rows2,
+ set_idx,
+ toList,
+ transpose,
+ zip2,
+ zipwith2 )
+import Linear.System ( solve_positive_definite )
+import Piecewise ( Piecewise(..), from_intervals )
import Polynomials.Orthogonal ( legendre )
-- | Dirichlet boundary conditions. Since u(a)=u(b)=0 are fixed,
two = fromInteger 2
+-- * Load vector
+
-- | Normalized integrals of orthogonal basis functions over
--- n[-1,1]. The test case below comes from Sage where the
+-- [-1,1]. The test case below comes from Sage where the
-- orthogonality of the polynomials' derivatives can easily be
-- tested.
--
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
-- 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.
(Arity l, Arity m, Arity n,
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
--
-- 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
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
+
+
+solution :: forall m n l a.
+ (Arity m, Arity n, Arity l,
+ Algebraic.C a, Eq a, RealField.C a, ToRational.C a, Show a)
+ => PDE a
+ -> Params m n (S l) a
+ -> Piecewise a
+solution pde params =
+ from_intervals $ map head $ toList $ solved_column
+ where
+ global_coeffs :: Col (S l) a
+ global_coeffs = coefficients pde params
+
+ ptr :: Mat m (S n) Int
+ ptr = pointer params
+
+ -- Each mesh element has an associated row in the pointer
+ -- matrix. Stick them together.
+ mesh_with_ptr_rows :: Col m (Interval a, Row (S n) Int)
+ mesh_with_ptr_rows = zip2 (mesh params) (rows2 ptr)
+
+ make_local_coeffs :: (Interval a, Row (S n) Int) -> Row (S n) a
+ make_local_coeffs (interval, ptr_row) =
+ construct lambda
+ where
+ lambda _ j = if (ptr_row !!! (0,j)) == zero
+ then zero
+ else global_coeffs !!! ((ptr_row !!! (0,j)) - 1, 0)
+
+ -- Create a column vector for each mesh element containing the global
+ -- coefficients corresponding to that element.
+ local_coeffs :: Col m (Row (S n) a)
+ local_coeffs = map2 make_local_coeffs mesh_with_ptr_rows
+
+ global_basis_functions :: Col (S n) (a -> a)
+ global_basis_functions =
+ construct lambda
+ where lambda i _ = big_N (toInteger i)
+
+ mesh_with_coeffs :: Col m (Interval a, Row (S n) a)
+ mesh_with_coeffs = zip2 (mesh params) local_coeffs
+
+ solved_column :: Col m (Interval a, (a -> a))
+ solved_column = map2 solve_piece $ mesh_with_coeffs
+
+ solve_piece :: (Interval a, Row (S n) a) -> (Interval a, (a -> a))
+ solve_piece (interval, coeffs_row) = (interval, f)
+ where
+ coeffs_col = transpose coeffs_row
+
+ f x = element_sum2 $ zipwith2 combine coeffs_col global_basis_functions
+ where
+ xi = (affine interval) x
+ combine ci ni = ci*(ni xi)
+
+
+
+
+solution' :: forall m n l a.
+ (Arity m, Arity n, Arity l,
+ Algebraic.C a, Eq a, RealField.C a, ToRational.C a, Show a)
+ => Col (S l) a
+ -> PDE a
+ -> Params m n (S l) a
+ -> Piecewise a
+solution' global_coeffs pde params =
+ from_intervals $ map head $ toList $ solved_column
+ where
+-- global_coeffs :: Col (S l) a
+-- global_coeffs = coefficients pde params
+
+ ptr :: Mat m (S n) Int
+ ptr = pointer params
+
+ -- Each mesh element has an associated row in the pointer
+ -- matrix. Stick them together.
+ mesh_with_ptr_rows :: Col m (Interval a, Row (S n) Int)
+ mesh_with_ptr_rows = zip2 (mesh params) (rows2 ptr)
+
+ make_local_coeffs :: (Interval a, Row (S n) Int) -> Row (S n) a
+ make_local_coeffs (interval, ptr_row) =
+ construct lambda
+ where
+ lambda _ j = if (ptr_row !!! (0,j)) == zero
+ then zero
+ else global_coeffs !!! ((ptr_row !!! (0,j)) - 1, 0)
+
+ -- Create a column vector for each mesh element containing the global
+ -- coefficients corresponding to that element.
+ local_coeffs :: Col m (Row (S n) a)
+ local_coeffs = map2 make_local_coeffs mesh_with_ptr_rows
+
+ global_basis_functions :: Col (S n) (a -> a)
+ global_basis_functions =
+ construct lambda
+ where lambda i _ = big_N (toInteger i)
+
+ mesh_with_coeffs :: Col m (Interval a, Row (S n) a)
+ mesh_with_coeffs = zip2 (mesh params) local_coeffs
+
+ solved_column :: Col m (Interval a, (a -> a))
+ solved_column = map2 solve_piece $ mesh_with_coeffs
+
+ solve_piece :: (Interval a, Row (S n) a) -> (Interval a, (a -> a))
+ solve_piece (interval, coeffs_row) = (interval, f)
+ where
+ coeffs_col = transpose coeffs_row
+
+ f x = element_sum2 $ zipwith2 combine coeffs_col global_basis_functions
+ where
+ xi = (affine interval) x
+ combine ci ni = ci*(ni xi)
+
projects / numerical-analysis.git / blobdiff