+
+
+solution :: 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
+ -> Piecewise a
+solution pde params =
+ from_intervals $ map head $ toList $ solved_column
+ where
+ global_coeffs :: Col (S l) a
+ global_coeffs = coefficients pde params
+
+ ptr :: Mat m (S n) Int
+ ptr = pointer params
+
+ -- Each mesh element has an associated row in the pointer
+ -- matrix. Stick them together.
+ mesh_with_ptr_rows :: Col m (Interval a, Row (S n) Int)
+ mesh_with_ptr_rows = zip2 (mesh params) (rows2 ptr)
+
+ make_local_coeffs :: (Interval a, Row (S n) Int) -> Row (S n) a
+ make_local_coeffs (_, ptr_row) =
+ construct lambda
+ where
+ lambda _ j = if (ptr_row !!! (0,j)) == zero
+ then zero
+ else global_coeffs !!! ((ptr_row !!! (0,j)) - 1, 0)
+
+ -- Create a column vector for each mesh element containing the global
+ -- coefficients corresponding to that element.
+ local_coeffs :: Col m (Row (S n) a)
+ local_coeffs = map2 make_local_coeffs mesh_with_ptr_rows
+
+ global_basis_functions :: Col (S n) (a -> a)
+ global_basis_functions =
+ construct lambda
+ where lambda i _ = big_N (toInteger i)
+
+ mesh_with_coeffs :: Col m (Interval a, Row (S n) a)
+ mesh_with_coeffs = zip2 (mesh params) local_coeffs
+
+ solved_column :: Col m (Interval a, (a -> a))
+ solved_column = map2 solve_piece $ mesh_with_coeffs
+
+ solve_piece :: (Interval a, Row (S n) a) -> (Interval a, (a -> a))
+ solve_piece (interval, coeffs_row) = (interval, g)
+ where
+ coeffs_col = transpose coeffs_row
+
+ g x = element_sum2 $ zipwith2 combine coeffs_col global_basis_functions
+ 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
+
+
+
+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