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 Integration.Gaussian ( gaussian )
import Linear.Matrix (
Col,
Mat(..),
(!!!),
construct,
+ ifoldl2,
nrows )
import Polynomials.Orthogonal ( legendre )
-- | Dirichlet boundary conditions. Since u(a)=u(b)=0 are fixed,
-- there's no additional information conveyed by this type.
-data Dirichlet = Dirichlet
+data Dirichlet a = Dirichlet { domain_dirichlet :: Interval a }
-- | Neumann boundary conditions. @alpha@ specifies A(a)u'(b) and
-- @beta@ specifies A(b)u'(b).
-data Neumann a = Neumann { alpha :: a, beta :: a }
+data Neumann a =
+ Neumann { domain_neumann :: Interval a,
+ alpha :: a,
+ beta :: a }
-- | Boundary conditions can be either Dirichlet or Neumann.
-type BoundaryConditions a = Either Dirichlet (Neumann a)
+type BoundaryConditions a = Either (Dirichlet a) (Neumann a)
type Interval a = (a,a)
c :: (a -> a),
-- | f(x)
f :: (a -> a),
- -- | The domain in R^1 as an interval
- domain :: Interval a,
+
+ -- | The boundary conditions. The domain also specifies the
+ -- boundary in R^1.
bdy :: BoundaryConditions a }
-- | Non-PDE parameters for the finite element method. The additional
-- type parameter @n@ should be a type-level representation of the
-- largest element in @max_degrees@. It needs to be known statically
--- for the dimensions of the pointer matrix.
-data Params m n a =
+-- for the dimensions of the pointer matrix. The parameter @l@ is
+-- the number of global basis functions. It's equal to the number of
+-- /internal/ mesh nodes (i.e. m-1), plus the sum of (p_i - 1) for
+-- each p_i in max_degrees.
+data Params m n l a =
Params {
-- | A partition of the domain.
mesh :: Col m (Interval a),
--
-- >>> import Data.Vector.Fixed ( N5, N6 )
-- >>> import Linear.Matrix ( Col5, fromList )
+-- >>> import Naturals ( N19 )
--
-- >>> let p = fromList [[3],[3],[5],[4],[5]] :: Col5 Int
-- >>> let mesh = undefined :: Col5 (Int,Int)
--- >>> let params = Params mesh p :: Params N5 N5 Int
+-- >>> let params = Params mesh p :: Params N5 N5 N19 Int
-- >>> let row1 = [0,1,5,6,0,0] :: [Int]
-- >>> let row2 = [1,2,7,8,0,0] :: [Int]
-- >>> let row3 = [2,3,9,10,11,12] :: [Int]
-- >>> pointer params == expected
-- True
--
-pointer :: (Arity m, Arity n) => Params m n a -> Mat m (S n) Int
+pointer :: (Arity m, Arity n, Arity l) => Params m n l a -> Mat m (S n) Int
pointer params =
construct lambda
where
--
-- Examples:
--
--- >>> let phi = affine (-6,9)
+-- >>> let phi = affine (-6,7)
-- >>> phi (-6)
-- -1.0
--- >>> phi (9)
+-- >>> phi 7
-- 1.0
--
affine :: Field.C a => (a,a) -> (a -> a)
affine (x1,x2) x = (fromInteger 2)*(x - x1)/(x2 - x1) - (fromInteger 1)
+-- | The inverse of 'affine'. It should send [-1,1] into [x1,x2].
+--
+-- Examples:
+--
+-- >>> let phi = affine_inv (-6,7)
+-- >>> phi (-1)
+-- -6.0
+-- >>> phi 1
+-- 7.0
+--
+affine_inv :: Field.C a => (a,a) -> (a -> a)
+affine_inv (x1,x2) x =
+ x*(x2 - x1)/two + (x1 + x2)/two
+ where
+ two = fromInteger 2
+
-- | Orthonormal basis functions over [-1,1]. The test case below
-- comes from Sage where the orthogonality of the polynomials'
--
-- >>> import qualified Algebra.Absolute as Absolute ( abs )
--
--- >>> let expected = 6.33910180790284
--- >>> let actual = big_N 3 1.5 :: Double
+-- >>> let expected = 2.99624907925257
+-- >>> let actual = big_N 4 1.5 :: Double
-- >>> Absolute.abs (actual - expected) < 1e-12
-- True
--
big_N :: forall a. (Algebraic.C a, RealField.C a) => Integer -> a -> a
big_N k x
| k < 0 = error "requested a negative basis function"
+ | k == 0 = (one - x) / (fromInteger 2)
+ | k == 1 = (one + x) / (fromInteger 2)
| otherwise =
- coeff * ( legendre (k+2) x - legendre k x )
+ coeff * ( legendre k x - legendre (k-2) x )
where
+ two = fromInteger 2
four = fromInteger 4
- six = fromInteger 6
- coeff = one / (sqrt (four*(fromInteger k) + six)) :: a
+ coeff = one / (sqrt (four*(fromInteger k) - two)) :: a
+
+
+-- | A matrix containing 'big_N' functions indexed by their
+-- element/number. Each row in the matrix represents a finite element
+-- (i.e. an interval in the mesh). Within row @i@, column @j@ contains
+-- the @j@th 'big_N' basis function.
+--
+-- Any given 'big_N' will probably wind up in this matrix multiple
+-- times; the basis functions don't change depending on the
+-- interval. Only the /number/ of basis functions does. Computing
+-- them this way allows us to easily construct a lookup \"table\" of
+-- the proper dimensions.
+--
+-- The second example below relies on the fact that @big_N 3@ and
+-- @big_N 6@ expand to Legendre polynomials (2,4) and (5,7)
+-- respectively and so should be orthogonal over [-1,1].
+--
+-- Examples:
+--
+-- >>> import Data.Vector.Fixed ( N5 )
+-- >>> import Integration.Gaussian ( gaussian )
+-- >>> import Linear.Matrix ( Col5, fromList )
+-- >>> import Naturals ( N19 )
+--
+-- >>> let p = fromList [[3],[3],[5],[4],[5]] :: Col5 Int
+-- >>> let mesh = undefined :: Col5 (Double,Double)
+-- >>> let params = Params mesh p :: Params N5 N5 N19 Double
+-- >>> let big_ns = big_N_matrix params
+-- >>> let n1 = big_ns !!! (1,0)
+-- >>> let n4 = big_ns !!! (4,0)
+-- >>> n1 1.5 == n4 1.5
+-- True
+-- >>> let n1 = big_ns !!! (1,3)
+-- >>> let n2 = big_ns !!! (2,4)
+-- >>> gaussian (\x -> (n1 x) * (n2 x)) < 1e-12
+-- True
+--
+big_N_matrix :: (Arity m, Arity n, Arity l, Algebraic.C a, RealField.C a)
+ => Params m n l a
+ -> Mat m (S n) (a -> a)
+big_N_matrix params =
+ construct lambda
+ where
+ maxdeg i = (max_degrees params) !!! (i,0)
+ lambda i j = if j > maxdeg i
+ then (const $ fromInteger 0)
+ else big_N (toInteger j)
+
+
+
+-- | Compute the global load vector F.
+--
+-- Examples:
+--
+-- >>> import Data.Vector.Fixed ( N3, N4 )
+-- >>> import Linear.Matrix ( Col4, frobenius_norm, fromList )
+-- >>> import Naturals ( N7 )
+--
+-- >>> let big_A = const (1::Double)
+-- >>> let c x = sin x
+-- >>> let f x = x*(sin x)
+-- >>> let bdy = Left (Dirichlet (0,1::Double))
+-- >>> let pde = PDE big_A c f bdy
+--
+-- >>> let i1 = (0.0,1/3)
+-- >>> let i2 = (1/3,2/3)
+-- >>> let i3 = (2/3,4/5)
+-- >>> let i4 = (4/5,1.0)
+-- >>> let mesh = fromList [[i1], [i2], [i3], [i4]] :: Col4 (Double,Double)
+-- >>> let pvec = fromList [[2],[3],[2],[1]] :: Col4 Int
+-- >>> let params = Params mesh pvec :: Params N4 N3 N7 Double
+--
+-- >>> let f1 = [0.0418]
+-- >>> let f2 = [0.0805]
+-- >>> let f3 = [0.1007]
+-- >>> let f4 = [-0.0045]
+-- >>> let f5 = [-0.0332]
+-- >>> let f6 = [-0.0054]
+-- >>> let f7 = [-0.0267]
+-- >>> let expected = fromList [f1,f2,f3,f4,f5,f6,f7] :: Col N7 Double
+-- >>> let actual = big_F pde params
+-- >>> frobenius_norm (actual - expected) < 1e-4
+-- True
+--
+big_F :: forall m n l a.
+ (Arity l, Arity m, Arity n,
+ Algebraic.C a, RealField.C a, ToRational.C a)
+ => PDE a
+ -> Params m n l a
+ -> Col l a
+big_F pde params =
+ ifoldl2 accum zero (big_N_matrix params)
+ where
+ accum :: Int -> Int -> Col l a -> (a -> a) -> Col l a
+ accum i j prev_F this_N =
+ prev_F + this_F
+ where
+ two = fromInteger 2
+ (x1,x2) = (mesh params) !!! (i,0)
+ q = affine_inv (x1,x2)
+ integrand x = ((f pde) (q x)) * (this_N x)
+
+ -- The pointer matrix numbers from 1 so subtract one here to
+ -- get the right index.
+ global_idx = ((pointer params) !!! (i,j)) - 1
+ this_F = construct lambda
+ lambda k _ = if k == global_idx
+ then (gaussian integrand)*(x2 - x1) / two
+ else fromInteger 0
projects / numerical-analysis.git / commitdiff