module FEM.R1
where
+import Algebra.Absolute ( abs )
import qualified Algebra.Algebraic as Algebraic ( C )
import qualified Algebra.Field as Field ( C )
import qualified Algebra.RealField as RealField ( C )
import qualified Algebra.ToRational as ToRational ( C )
import Data.Vector.Fixed ( Arity, S )
-import NumericPrelude
+import NumericPrelude hiding ( abs )
import qualified Prelude as P
import Integration.Gaussian ( gaussian )
Row,
(!!!),
construct,
+ dot,
+ element_sum2,
fromList,
ifoldl2,
+ map2,
nrows,
- set_idx )
+ rows2,
+ set_idx,
+ toList,
+ transpose,
+ zip2,
+ zipwith2 )
import Linear.System ( solve_positive_definite )
+import Piecewise ( Piecewise(..), evaluate', from_intervals )
import Polynomials.Orthogonal ( legendre )
-- | Dirichlet boundary conditions. Since u(a)=u(b)=0 are fixed,
-- * 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.
--
--
-- 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
+-- >>> import Linear.Matrix ( Col7, frobenius_norm )
+-- >>> import FEM.R1.Example1 ( pde', params' )
--
-- >>> let f1 = [0.0418]
-- >>> let f2 = [0.0805]
-- >>> 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
+-- >>> let expected = fromList [f1,f2,f3,f4,f5,f6,f7] :: Col7 Double
+-- >>> let actual = big_F pde' params'
-- >>> frobenius_norm (actual - expected) < 1e-4
-- True
--
-> 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))
+ ifoldl2 accum cur_K (big_N's_matrix :: Mat (S n) (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 =
(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.
+ -- The pointer matrix numbers from 1 so subtract one below to
+ -- get the right index. The indices i,j have upper bounds
+ -- dependent on the element k. Since we statically create the
+ -- matrix of basis function derivatives, we have to check here
+ -- whether or not i,j exceed the max index.
row_idx = ((pointer params) !!! (k,i)) - 1
col_idx = ((pointer params) !!! (k,j)) - 1
integral = (two/(x2 - x1))* (gaussian integrand)
--
-- 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
+-- >>> import Linear.Matrix ( Mat7, frobenius_norm )
+-- >>> import FEM.R1.Example1 ( pde', params' )
--
-- >>> let k1 = [6, -3, 0, 0, 0, 0, 0] :: [Double]
-- >>> let k2 = [-3, 10.5, -7.5, 0, 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
+-- >>> let expected = fromList [k1,k2,k3,k4,k5,k6,k7] :: Mat7 Double
+-- >>> let actual = big_K pde' params'
-- >>> frobenius_norm (actual - expected) < 1e-10
-- True
--
-> 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))
+ ifoldl2 accum cur_M (big_Ns_matrix :: Mat (S n) (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 =
--
-- 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
+-- >>> import Linear.Matrix ( Mat7, frobenius_norm )
+-- >>> import FEM.R1.Example1 ( pde', params' )
--
-- >>> 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 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
+-- >>> let expected = fromList [m1,m2,m3,m4,m5,m6,m7] :: Mat7 Double
+-- >>> let actual = big_M pde' params'
-- >>> frobenius_norm (actual - expected) < 1e-3
-- True
--
--
-- 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
+-- >>> import Linear.Matrix ( Col7, frobenius_norm )
+-- >>> import FEM.R1.Example1 ( pde', params' )
--
-- >>> let c1 = [0.02366220347687] :: [Double]
-- >>> let c2 = [0.03431630082636] :: [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
+-- >>> let actual = coefficients pde' params'
-- >>> frobenius_norm (actual - expected) < 1e-8
-- True
--
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)
+ => 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 (_, 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, g)
+ where
+ coeffs_col = transpose coeffs_row
+
+ g x = element_sum2 $ zipwith2 combine coeffs_col global_basis_functions
+ where
+ xi = (affine interval) x
+ combine ci ni = ci*(ni xi)
+
+
+-- energy_true :: (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
+-- -> (a -> a) -- ^ True solution @u@
+-- -> (a -> a) -- ^ Derivative of true solution @u'@
+-- -> a
+-- energy_true pde params u u' =
+-- case (bdy pde) of
+-- Left (Dirichlet (x1,x2)) ->
+-- sqrt $ bilinear_form u u' u u'
+-- where
+-- two = fromInteger 2
+-- q = affine_inv (x1,x2)
+-- bilinear_form w w' v v' = (x2 - x1)*(gaussian integrand)/two
+-- where
+-- integrand x = ((big_A pde) (q x))*(w' (q x))*(v' (q x))
+-- + ((c pde) (q x))*(w (q x))*(v (q x))
+
+-- _ -> error "Neumann BCs not implemented."
+
+
+energy_fem :: (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
+ -> a
+energy_fem pde params =
+ (coefficients pde params) `dot` (big_F pde params)
+
+
+relative_error :: 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
+ -> a -- ^ The energy norm of the true solution @u@
+ -> a
+relative_error pde params energy_true =
+ cent * sqrt(energy_true - (energy_fem pde params)/energy_true)
+ where
+ cent = fromInteger 100
+
+
+
+relative_error_pointwise :: 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
+ -> (a -> a) -- ^ The true solution @u@
+ -> a -- ^ The point @x@ at which to compute the error.
+ -> a
+relative_error_pointwise pde params u x =
+ cent * ( u_exact - u_fem ) / u_exact
+ where
+ u_exact = abs $ u x
+ u_fem = evaluate' (solution pde params) x
+ cent = fromInteger 100
projects / numerical-analysis.git / blobdiff