X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FFEM%2FR1.hs;h=7b89ffabdf41513197e77d93f6952fd09dc3706d;hb=84bf90406341348a41c7cc57c3d365238690a930;hp=f99b3dcace7be32ee7df3ba07a488dc012285ae2;hpb=d0056767bcf231e1d5cce78a85f6be0688770cc3;p=numerical-analysis.git

diff --git a/src/FEM/R1.hs b/src/FEM/R1.hs
index f99b3dc..7b89ffa 100644
--- a/src/FEM/R1.hs
+++ b/src/FEM/R1.hs
@@ -35,7 +35,7 @@ import qualified Algebra.RealField as RealField ( C )
 import qualified Algebra.ToRational as ToRational ( C )
 import Data.Vector.Fixed ( Arity, S )
 import NumericPrelude hiding ( abs )
-import qualified Prelude as P
+import Prelude ()
 
 import Integration.Gaussian ( gaussian )
 import Linear.Matrix (
@@ -190,11 +190,11 @@ affine (x1,x2) x = (fromInteger 2)*(x - x1)/(x2 - x1) - (fromInteger 1)
 --   >>> phi 1
 --   7.0
 --
-affine_inv :: Field.C a => (a,a) -> (a -> a)
+affine_inv :: forall a. Field.C a => (a,a) -> (a -> a)
 affine_inv (x1,x2) x =
   x*(x2 - x1)/two + (x1 + x2)/two
   where
-    two = fromInteger 2
+    two = fromInteger 2 :: a
 
 
 -- * Load vector
@@ -221,9 +221,9 @@ big_N k x
   | otherwise =
       coeff * ( legendre k x - legendre (k-2) x )
       where
-        two = fromInteger 2
-        four = fromInteger 4
-        coeff = one / (sqrt (four*(fromInteger k) - two)) :: a
+        two = fromInteger 2 :: a
+        four = fromInteger 4 :: a
+        coeff = one / (sqrt (four*(fromInteger k) - two))
 
 
 -- | A matrix containing 'big_N' functions indexed by their
@@ -311,7 +311,7 @@ big_F pde params =
     accum i j prev_F this_N =
       prev_F + this_F
       where
-        two = fromInteger 2
+        two = fromInteger 2 :: a
         (x1,x2) = (mesh params) !!! (i,0)
         q = affine_inv (x1,x2)
         integrand x = ((f pde) (q x)) * (this_N x)
@@ -347,8 +347,8 @@ big_N' k x
   | k == 1 = one / (fromInteger 2)
   | otherwise = coeff * ( legendre k x )
       where
-        two = fromInteger 2
-        coeff = sqrt ((two*(fromInteger k) + one) / two) :: a
+        two = fromInteger 2 :: a
+        coeff = sqrt ((two*(fromInteger k) + one) / two)
 
 
 -- | The matrix of (N_i' * N_j') functions used in the integrand of
@@ -378,7 +378,7 @@ big_K_elem pde params _ k cur_K _ =
     accum i j prev_K these_N's =
       prev_K + this_K
       where
-        two = fromInteger 2
+        two = fromInteger 2 :: a
         (x1,x2) = (mesh params) !!! (k,0)
         q = affine_inv (x1,x2)
         integrand x = ((big_A pde) (q x)) * (these_N's x)
@@ -447,7 +447,7 @@ big_M_elem pde params _ k cur_M _ =
     accum i j prev_M these_Ns =
       prev_M + this_M
       where
-        two = fromInteger 2
+        two = fromInteger 2 :: a
         (x1,x2) = (mesh params) !!! (k,0)
         q = affine_inv (x1,x2)
         integrand x = ((c pde) (q x)) * (these_Ns x)
@@ -583,33 +583,11 @@ solution pde params =
             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
+           => PDE a
+           -> Params m n (S l) a
+           -> a
 energy_fem pde params =
   (coefficients pde params) `dot` (big_F pde params)
 
@@ -624,7 +602,7 @@ relative_error :: forall m n l a.
 relative_error pde params energy_true =
   cent * sqrt(energy_true - (energy_fem pde params)/energy_true)
   where
-    cent = fromInteger 100
+    cent = fromInteger 100 :: a
 
 
 
@@ -637,8 +615,7 @@ relative_error_pointwise :: forall m n l a.
                          -> 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
+  cent * ( abs $ (u x) - u_fem ) / ( abs $ u x )
   where
-    u_exact = abs $ u x
     u_fem = evaluate' (solution pde params) x
-    cent = fromInteger 100
+    cent = fromInteger 100 :: a