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 qualified Prelude as P
+import NumericPrelude hiding ( abs )
+import Prelude ()
import Integration.Gaussian ( gaussian )
import Linear.Matrix (
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,
-- >>> 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
-> 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 =
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)
- -- 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)
-> 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 =
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)
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_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 :: a
+
+
+
+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 * ( abs $ (u x) - u_fem ) / ( abs $ u x )
+ where
+ u_fem = evaluate' (solution pde params) x
+ cent = fromInteger 100 :: a
projects / numerical-analysis.git / blobdiff