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 (
-- >>> 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
| 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
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)
| 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
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)
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)
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)
relative_error pde params energy_true =
cent * sqrt(energy_true - (energy_fem pde params)/energy_true)
where
- cent = fromInteger 100
+ cent = fromInteger 100 :: 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