X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FFEM%2FR1.hs;h=28317a6faecabaf0a079215a192f961a7d5e5be3;hb=45aa63a662556bc9ed0f6018f9f3f580586f38a9;hp=217abd3009acaac271e2a17fe99eccd6e5df9050;hpb=25706b7e253af82d8ca9128eb23772c47ee6c0c3;p=numerical-analysis.git

diff --git a/src/FEM/R1.hs b/src/FEM/R1.hs
index 217abd3..28317a6 100644
--- a/src/FEM/R1.hs
+++ b/src/FEM/R1.hs
@@ -28,12 +28,13 @@
 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 )
@@ -43,6 +44,7 @@ import Linear.Matrix (
   Row,
   (!!!),
   construct,
+  dot,
   element_sum2,
   fromList,
   ifoldl2,
@@ -55,7 +57,7 @@ import Linear.Matrix (
   zip2,
   zipwith2 )
 import Linear.System ( solve_positive_definite )
-import Piecewise ( Piecewise(..), from_intervals )
+import Piecewise ( Piecewise(..), evaluate', from_intervals )
 import Polynomials.Orthogonal ( legendre )
 
 -- | Dirichlet boundary conditions. Since u(a)=u(b)=0 are fixed,
@@ -281,22 +283,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]
@@ -305,8 +293,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
 --
@@ -384,7 +372,7 @@ big_K_elem :: forall m n l a b.
       -> 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 =
@@ -394,8 +382,11 @@ big_K_elem pde params _ k cur_K _ =
         (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)
@@ -410,22 +401,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]
@@ -434,8 +411,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
 --
@@ -464,7 +441,7 @@ big_M_elem :: forall m n l a b.
       -> 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 =
@@ -489,22 +466,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]
@@ -514,8 +477,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
 --
@@ -536,22 +499,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]
@@ -560,8 +509,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
 --
@@ -580,7 +530,7 @@ coefficients 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)
+            Algebraic.C a, Eq a, RealField.C a, ToRational.C a)
          => PDE a
          -> Params m n (S l) a
          -> Piecewise a
@@ -599,7 +549,7 @@ solution pde params =
     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) =
+    make_local_coeffs (_, ptr_row) =
       construct lambda
         where
           lambda _ j = if (ptr_row !!! (0,j)) == zero
@@ -623,70 +573,50 @@ solution pde params =
     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)
+    solve_piece (interval, coeffs_row) = (interval, g)
       where
         coeffs_col = transpose coeffs_row
 
-        f x = element_sum2 $ zipwith2 combine coeffs_col global_basis_functions
+        g 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
+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
---    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)
+    cent = fromInteger 100
 
-    -- 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)
 
+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