X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FFEM%2FR1.hs;h=ae7548a86d75d42311691c8bdcfa780384b3938c;hb=3a58695c57f57c9adbc0fb2ca26e06e7da7bfd63;hp=922469f730dff4af6ed9d8387928c60335acd223;hpb=1e988aee439f5352c25a4de03d7d18a58ae1d403;p=numerical-analysis.git diff --git a/src/FEM/R1.hs b/src/FEM/R1.hs index 922469f..ae7548a 100644 --- a/src/FEM/R1.hs +++ b/src/FEM/R1.hs @@ -43,11 +43,19 @@ import Linear.Matrix ( Row, (!!!), construct, + element_sum2, fromList, ifoldl2, + map2, nrows, - set_idx ) + 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, @@ -190,7 +198,7 @@ affine_inv (x1,x2) x = -- * 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. -- @@ -273,22 +281,8 @@ big_Ns_matrix = -- -- 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] @@ -297,8 +291,8 @@ big_Ns_matrix = -- >>> 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 -- @@ -402,22 +396,8 @@ big_K_elem pde params _ k cur_K _ = -- -- 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] @@ -426,8 +406,8 @@ big_K_elem pde params _ k cur_K _ = -- >>> 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 -- @@ -481,22 +461,8 @@ big_M_elem pde params _ k cur_M _ = -- -- 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] @@ -506,8 +472,8 @@ big_M_elem pde params _ k cur_M _ = -- >>> 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 -- @@ -528,22 +494,8 @@ big_M pde params = -- -- 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] @@ -552,8 +504,9 @@ big_M pde params = -- >>> 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 -- @@ -568,3 +521,58 @@ coefficients pde params = 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 (_, 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)