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 )
Row,
(!!!),
construct,
+ dot,
element_sum2,
fromList,
ifoldl2,
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,
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
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