Add relative error functions to FEM.R1.

diff --git a/src/FEM/R1.hs b/src/FEM/R1.hs
index 3e6869557d66f5601c769ff3863448a666729d60..f99b3dcace7be32ee7df3ba07a488dc012285ae2 100644 (file)
--- 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,
@@ -528,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
@@ -579,3 +581,64 @@ solution pde params =
           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