clifford 0.1.0.1 → 0.1.0.3
raw patch · 6 files changed
+116/−100 lines, 6 files
Files
- clifford.cabal +3/−3
- src/Numeric/Clifford/Blade.lhs +37/−31
- src/Numeric/Clifford/ClassicalMechanics.lhs +15/−15
- src/Numeric/Clifford/Multivector.lhs +34/−29
- src/Numeric/Clifford/NumericIntegration.lhs +21/−21
- src/Numeric/Clifford/NumericIntegration/DefaultIntegrators.hs +6/−1
clifford.cabal view
@@ -10,7 +10,7 @@ -- PVP summary: +-+------- breaking API changes -- | | +----- non-breaking API additions -- | | | +--- code changes with no API change-version: 0.1.0.1+version: 0.1.0.3 -- A short (one-line) description of the package. synopsis: A Clifford algebra library@@ -48,6 +48,7 @@ -- Constraint on the version of Cabal needed to build this package. cabal-version: >=1.10 + library -- Modules exported by the library. exposed-modules: Numeric.Clifford.Blade, Numeric.Clifford.Multivector, Numeric.Clifford.NumericIntegration, Numeric.Clifford.NumericIntegration.DefaultIntegrators, Numeric.Clifford.ClassicalMechanics@@ -56,8 +57,7 @@ -- other-modules: -- LANGUAGE extensions used by modules in this package.- -- other-extensions:- + other-extensions: NoImplicitPrelude, ScopedTypeVariables, GADTs, DataKinds,KindSignatures,UnicodeSyntax,FlexibleContexts, RankNTypes,TemplateHaskell,NoMonomorphismRestriction,MultiParamTypeClasses,FlexibleInstances, TypeOperators -- Other library packages from which modules are imported. build-depends: base >=4.6 && <4.9, numeric-prelude >= 0.4.0.1 && < 0.5.0, permutation >= 0.4.1 && < 0.5, data-ordlist >= 0.4.5 && < 0.5, converge >= 0.1.0.1 && < 0.2, lens >= 4.0.3 && < 4.1,
src/Numeric/Clifford/Blade.lhs view
@@ -49,8 +49,8 @@ import Test.QuickCheck import Control.Lens hiding (indices) import Data.DeriveTH-import GHC.TypeLits hiding (isEven)-import GHC.Real (fromIntegral)+import GHC.TypeLits hiding (isEven, isOdd)+import GHC.Real (fromIntegral, toInteger) import Algebra.Field import Debug.Trace --trace _ a = a@@ -63,19 +63,19 @@ \texttt{bScale} is the amplitude of the blade. \texttt{bIndices} are the indices for the basis. \begin{code} -data Blade (n :: Nat) f where- Blade :: forall n f . (SingI n, Algebra.Field.C f) => {_scale :: f, _indices :: [Natural]} -> Blade n f+data Blade (p :: Nat) (q :: Nat) f where+ Blade :: forall p q f . (SingI p, SingI q, Algebra.Field.C f) => {_scale :: f, _indices :: [Natural]} -> Blade p q f -- makeLenses ''Blade-scale :: Lens' (Blade n f) f+scale :: Lens' (Blade p q f) f scale = lens _scale (\blade v -> blade {_scale = v})-indices :: Lens' (Blade n f) [Natural]+indices :: Lens' (Blade p q f) [Natural] indices = lens _indices (\blade v -> blade {_indices = v})-dimension :: forall (n::Nat) f. SingI n => Blade n f -> Natural-dimension _ = toNatural ((GHC.Real.fromIntegral $ fromSing (sing :: Sing n))::Word)+dimension :: forall (p::Nat) (q::Nat) f. (SingI p, SingI q) => Blade p q f -> (Natural,Natural)+dimension _ = (toNatural ((GHC.Real.fromIntegral $ fromSing (sing :: Sing p))::Word),toNatural((GHC.Real.fromIntegral $ fromSing (sing :: Sing q))::Word)) bScale b = b^.scale bIndices b = b^.indices-instance(Show f) => Show (Blade n f) where+instance(Show f) => Show (Blade p q f) where --TODO: Do this with HaTeX show (Blade scale indices) = pref ++ if null indices then "" else basis where pref = show scale@@ -84,7 +84,7 @@ vecced index = "\\vec{e_{" ++ show index ++ "}}" -instance (Algebra.Additive.C f, Eq f) => Eq (Blade n f) where+instance (Algebra.Additive.C f, Eq f) => Eq (Blade p q f) where a == b = aScale == bScale && aIndices == bIndices where (Blade aScale aIndices) = bladeNormalForm a (Blade bScale bIndices) = bladeNormalForm b@@ -93,10 +93,10 @@ For example, a scalar could be constructed like so: \texttt{Blade s []} \begin{code}-scalarBlade :: (Algebra.Field.C f, SingI n) => f -> Blade n f+scalarBlade :: (Algebra.Field.C f, SingI p, SingI q) => f -> Blade p q f scalarBlade d = Blade d [] -zeroBlade :: (Algebra.Field.C f, SingI n) => Blade n f+zeroBlade :: (Algebra.Field.C f, SingI p, SingI q) => Blade p q f zeroBlade = scalarBlade Algebra.Additive.zero bladeNonZero b = b^.scale /= Algebra.Additive.zero@@ -117,32 +117,38 @@ \begin{code}-bladeNormalForm :: forall (n::Nat) f. Blade n f -> Blade n f+bladeNormalForm :: forall (p::Nat) (q::Nat) f. Blade p q f -> Blade p q f bladeNormalForm (Blade scale indices) = result where- result = if (any (\i -> (GHC.Real.fromIntegral i) >n') indices) then trace "Blade contains vector with i > d" zeroBlade else Blade scale' uniqueSorted- n' = fromSing (sing :: Sing n)+ result = if (any (\i -> (GHC.Real.toInteger i) > d) indices) then trace "Blade contains vector with i > d" zeroBlade else Blade scale' uniqueSorted+ p' = (fromSing (sing :: Sing p)) :: Integer+ q' = (fromSing (sing :: Sing q)) :: Integer+ d = p' + q' numOfIndices = length indices (sorted, perm) = Data.Permute.sort numOfIndices indices- scale' = if isEven perm then scale else negate scale- uniqueSorted = removeDupPairs sorted+ scale' = if (isEven perm) /= (negated) then scale else negate scale+ (uniqueSorted,negated) = removeDupPairs [] sorted False where- removeDupPairs [] = []- removeDupPairs [x] = [x]- removeDupPairs (x:y:rest) | x == y = removeDupPairs rest- | otherwise = x : removeDupPairs (y:rest)+ removeDupPairs :: [Natural] -> [Natural] -> Bool -> ([Natural],Bool)+ removeDupPairs accum [] negated = (accum,negated)+ removeDupPairs accum [x] negated = (accum++[x],negated)+ removeDupPairs accum (x:y:rest) negated | x == y = + if GHC.Real.toInteger x > p' + then removeDupPairs accum rest (not negated)+ else removeDupPairs accum rest negated+ | otherwise = removeDupPairs (accum++[x]) (y:rest) negated \end{code} What is the grade of a blade? Just the number of indices. \begin{code}-grade :: Blade n f -> Integer-grade = toInteger . length . bIndices +grade :: Blade p q f -> Integer+grade = GHC.Real.toInteger . length . bIndices -bladeIsOfGrade :: Blade n f -> Integer -> Bool+bladeIsOfGrade :: Blade p q f -> Integer -> Bool blade `bladeIsOfGrade` k = grade blade == k -bladeGetGrade :: Integer -> Blade n f -> Blade n f+bladeGetGrade :: Integer -> Blade p q f -> Blade p q f bladeGetGrade k blade@(Blade _ _) = if blade `bladeIsOfGrade` k then blade else zeroBlade \end{code}@@ -152,9 +158,9 @@ First up for operations: Blade multiplication. This is no more than assembling orthogonal vectors into k-vectors. \begin{code}-bladeMul :: Blade n f -> Blade n f-> Blade n f-bladeMul x@(Blade _ _) y@(Blade _ _)= bladeNormalForm $ Blade (bScale x Algebra.Ring.* bScale y) (bIndices x ++ bIndices y)-multiplyBladeList :: (SingI n, Algebra.Field.C f) => [Blade n f] -> Blade n f+bladeMul :: Blade p q f -> Blade p q f-> Blade p q f+bladeMul x@(Blade _ _) y@(Blade _ _)= bladeNormalForm $ Blade (bScale x Algebra.Ring.* bScale y) (bIndices x ++ bIndices y) -- HAVE TO MAKE Q INDICES SQUARE NEGATIVE+multiplyBladeList :: (SingI p, SingI q, Algebra.Field.C f) => [Blade p q f] -> Blade p q f multiplyBladeList [] = error "Empty blade list!" multiplyBladeList (a:[]) = a multiplyBladeList a = bladeNormalForm $ Blade scale indices where@@ -167,7 +173,7 @@ Next up: The outer (wedge) product, denoted by $∧$ :3 \begin{code}-bWedge :: Blade n f -> Blade n f -> Blade n f+bWedge :: Blade p q f -> Blade p q f -> Blade p q f bWedge x y = bladeNormalForm $ bladeGetGrade k xy where k = grade x + grade y@@ -179,7 +185,7 @@ \begin{code}-bDot ::Blade n f -> Blade n f -> Blade n f+bDot ::Blade p q f -> Blade p q f -> Blade p q f bDot x y = bladeNormalForm $ bladeGetGrade k xy where k = Algebra.Absolute.abs $ grade x - grade y@@ -194,7 +200,7 @@ Now for linear combinations of (possibly different basis) blades. To start with, let's order basis blades: \begin{code}-instance (Algebra.Additive.C f, Ord f) => Ord (Blade n f) where+instance (Algebra.Additive.C f, Ord f) => Ord (Blade p q f) where compare a b | bIndices a == bIndices b = compare (bScale a) (bScale b) | otherwise = compare (bIndices a) (bIndices b)
src/Numeric/Clifford/ClassicalMechanics.lhs view
@@ -16,7 +16,7 @@ \begin{code} {-# LANGUAGE NoImplicitPrelude, FlexibleContexts, RankNTypes, ScopedTypeVariables, DeriveDataTypeable #-} {-# LANGUAGE NoMonomorphismRestriction, UnicodeSyntax, GADTs, KindSignatures, DataKinds #-}-{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FlexibleInstances, TypeOperators #-} {-# LANGUAGE TemplateHaskell #-} {-# LANGUAGE MultiParamTypeClasses #-} \end{code}@@ -59,9 +59,9 @@ import Debug.Trace --trace _ a = a -data EnergyMethod (d::Nat) f = Hamiltonian{ _dqs :: [DynamicSystem d f -> Multivector d f], _dps :: [DynamicSystem d f -> Multivector d f]}+data EnergyMethod (p::Nat) (q::Nat) f = Hamiltonian{ _dqs :: [DynamicSystem p q f -> Multivector p q f], _dps :: [DynamicSystem p q f -> Multivector p q f]} -data DynamicSystem (d::Nat) f = DynamicSystem {_time :: f, coordinates :: [Multivector d f], _momenta :: [Multivector d f], _energyFunction :: EnergyMethod d f, _projector :: DynamicSystem d f -> DynamicSystem d f}+data DynamicSystem (p::Nat) (q::Nat) f = DynamicSystem {_time :: f, coordinates :: [Multivector p q f], _momenta :: [Multivector p q f], _energyFunction :: EnergyMethod p q f, _projector :: DynamicSystem p q f -> DynamicSystem p q f} makeLenses ''EnergyMethod makeLenses ''DynamicSystem@@ -80,14 +80,14 @@ Now to make a physical object. \begin{code}-data ReferenceFrame (d::Nat) t = ReferenceFrame {basisVectors :: [Multivector d t]}-psuedoScalar' :: forall f (d::Nat). (Ord f, Algebra.Field.C f, SingI d) => ReferenceFrame d f -> Multivector d f+data ReferenceFrame (p::Nat) (q::Nat) t = ReferenceFrame {basisVectors :: [Multivector p q t]}+psuedoScalar' :: forall f (p::Nat) (q::Nat). (Ord f, Algebra.Field.C f, SingI p, SingI q) => ReferenceFrame p q f -> Multivector p q f psuedoScalar' = multiplyList . basisVectors-psuedoScalar :: forall (d::Nat) f. (Ord f, Algebra.Field.C f, SingI d) => Multivector d f-psuedoScalar = one `e` [1..(toNatural ((fromIntegral $ fromSing (sing :: Sing d))::Word))] ++ a `cross` b = (negate $ one)`e`[1,2,3] * (a ∧ b)-data PhysicalVector (d::Nat) t = PhysicalVector {dimension :: Natural, r :: Multivector d t, referenceFrame :: ReferenceFrame d t}+data PhysicalVector (p::Nat) (q::Nat) t = PhysicalVector {dimension :: Natural, r :: Multivector p q t, referenceFrame :: ReferenceFrame p q t} {-squishToDimension (PhysicalVector d (BladeSum terms) f) = PhysicalVector d r' f where r' = BladeSum terms' where terms' = terms & filter (\(Blade _ ind) -> all (\k -> k <= d) ind)@@ -95,14 +95,14 @@ r' = BladeSum terms' where terms' = terms & filter (\(Blade _ ind) -> all (\k -> k <= d) ind)-} -data RigidBody (d::Nat) f where- RigidBody:: (Algebra.Field.C f, Algebra.Module.C f (Multivector d f)) => {position :: PhysicalVector d f,- _momentum :: PhysicalVector d f,+data RigidBody (p::Nat) (q::Nat) f where+ RigidBody:: (Algebra.Field.C f, Algebra.Module.C f (Multivector p q f)) => {position :: PhysicalVector p q f,+ _momentum :: PhysicalVector p q f, _mass :: f,- _attitude :: PhysicalVector d f,- _angularMomentum :: PhysicalVector d f,- _inertia :: PhysicalVector d f- } -> RigidBody d f+ _attitude :: PhysicalVector p q f,+ _angularMomentum :: PhysicalVector p q f,+ _inertia :: PhysicalVector p q f+ } -> RigidBody p q f --makeLenses ''RigidBody doesn't actually work {- Things to do:
src/Numeric/Clifford/Multivector.lhs view
@@ -21,7 +21,7 @@ {-# LANGUAGE NoImplicitPrelude, FlexibleContexts, RankNTypes, ScopedTypeVariables, DeriveDataTypeable #-} {-# LANGUAGE NoMonomorphismRestriction, UnicodeSyntax, GADTs#-} {-# LANGUAGE FlexibleInstances, StandaloneDeriving, KindSignatures, DataKinds #-}-{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE TemplateHaskell, TypeOperators #-} {-# LANGUAGE MultiParamTypeClasses #-} \end{code} %if False@@ -83,17 +83,17 @@ A multivector is nothing but a linear combination of basis blades. \begin{code}-data Multivector (n::Nat) f where- BladeSum :: forall n f . (SingI n, Algebra.Field.C f, Ord f) => { _terms :: [Blade n f]} -> Multivector n f+data Multivector (p::Nat) (q::Nat) f where+ BladeSum :: forall p q f . (Ord f, Algebra.Field.C f, SingI p, SingI q) => { _terms :: [Blade p q f]} -> Multivector p q f -deriving instance Eq (Multivector n f)-deriving instance Ord (Multivector n f)-deriving instance (Show f) => Show (Multivector n f)+deriving instance Eq (Multivector p q f)+deriving instance Ord (Multivector p q f)+deriving instance (Show f) => Show (Multivector p q f) -dimension :: forall (n::Nat) f. SingI n => Multivector n f -> Natural-dimension _ = toNatural ((fromIntegral $ fromSing (sing :: Sing n))::Word)+dimension :: forall (p::Nat) (q::Nat) f. (SingI p, SingI q) => Multivector p q f -> (Natural,Natural)+dimension _ = (toNatural ((fromIntegral $ fromSing (sing :: Sing p))::Word),toNatural ((fromIntegral $ fromSing (sing :: Sing q))::Word)) -terms :: Lens' (Multivector n f) [Blade n f]+terms :: Lens' (Multivector p q f) [Blade p q f] terms = lens _terms (\bladeSum v -> bladeSum {_terms = v}) mvNormalForm (BladeSum terms) = BladeSum $ if null resultant then [scalarBlade Algebra.Additive.zero] else resultant where@@ -102,7 +102,7 @@ addLikeTerms' = sumLikeTerms . groupLikeTerms -groupLikeTerms ::Eq f => [Blade n f] -> [[Blade n f]]+groupLikeTerms ::Eq f => [Blade p q f] -> [[Blade p q f]] groupLikeTerms = groupBy (\a b -> a^.indices == b^.indices) compensatedSum' :: (Algebra.Additive.C f) => [f] -> f@@ -135,7 +135,7 @@ -- ei = multiplyAdd eim1 ai pii -multiplyOutBlades :: (SingI n, Algebra.Ring.C a) => [Blade n a] -> [Blade n a] -> [Blade n a]+multiplyOutBlades :: (SingI p, SingI q, Algebra.Ring.C a) => [Blade p q a] -> [Blade p q a] -> [Blade p q a] multiplyOutBlades x y = [bladeMul l r | l <-x, r <- y] --multiplyList :: Algebra.Ring.C t => [Multivector t] -> Multivector t@@ -150,29 +150,29 @@ sumList xs = mvNormalForm $ BladeSum $ concat $ map mvTerms xs -sumLikeTerms :: (Algebra.Field.C f, SingI n) => [[Blade n f]] -> [Blade n f]+sumLikeTerms :: (Algebra.Field.C f, SingI p, SingI q) => [[Blade p q f]] -> [Blade p q f] sumLikeTerms blades = map (\sameIxs -> map bScale sameIxs & compensatedSum' & (\result -> Blade result ((head sameIxs) & bIndices))) blades -instance (Algebra.Field.C f, SingI n, Ord f) => Data.Monoid.Monoid (Sum (Multivector n f)) where+instance (Algebra.Field.C f, SingI p, SingI q, Ord f) => Data.Monoid.Monoid (Sum (Multivector p q f)) where mempty = Data.Monoid.Sum Algebra.Additive.zero mappend (Data.Monoid.Sum x) (Data.Monoid.Sum y)= Data.Monoid.Sum (x + y) mconcat = Data.Monoid.Sum . sumList . map getSum -instance (Algebra.Field.C f, SingI n, Ord f) => Data.Monoid.Monoid (Product (Multivector n f)) where+instance (Algebra.Field.C f, SingI p, SingI q, Ord f) => Data.Monoid.Monoid (Product (Multivector p q f)) where mempty = Product one mappend (Product x) (Product y) = Product (x * y) mconcat = Product . foldl (*) one . map getProduct --Constructs a multivector from a scaled blade.-e :: (Algebra.Field.C f, Ord f, SingI n) => f -> [Natural] -> Multivector n f+e :: (Algebra.Field.C f, Ord f, SingI p, SingI q) => f -> [Natural] -> Multivector p q f s `e` indices = mvNormalForm $ BladeSum [Blade s indices] scalar s = s `e` [] -instance (Control.DeepSeq.NFData f) => Control.DeepSeq.NFData (Multivector n f)-instance (Control.DeepSeq.NFData f) => Control.DeepSeq.NFData (Blade n f)-instance (Algebra.Field.C f, Ord f, SingI n) => Algebra.Additive.C (Multivector n f) where+instance (Control.DeepSeq.NFData f) => Control.DeepSeq.NFData (Multivector p q f)+instance (Control.DeepSeq.NFData f) => Control.DeepSeq.NFData (Blade p q f)+instance (Algebra.Field.C f, Ord f, SingI p, SingI q) => Algebra.Additive.C (Multivector p q f) where a + b = mvNormalForm $ BladeSum (mvTerms a ++ mvTerms b) a - b = mvNormalForm $ BladeSum (mvTerms a ++ map bladeNegate (mvTerms b)) zero = BladeSum [scalarBlade Algebra.Additive.zero]@@ -184,7 +184,7 @@ \begin{code} -instance (Algebra.Field.C f, Ord f, SingI n) => Algebra.Ring.C (Multivector n f) where+instance (Algebra.Field.C f, Ord f, SingI p, SingI q) => Algebra.Ring.C (Multivector p q f) where BladeSum [Blade s []] * b = BladeSum $ map (bladeScaleLeft s) $ mvTerms b a * BladeSum [Blade s []] = BladeSum $ map (bladeScaleRight s) $ mvTerms a a * b = mvNormalForm $ BladeSum [bladeMul x y | x <- mvTerms a, y <- mvTerms b]@@ -199,6 +199,11 @@ two = fromInteger 2 mul = (Algebra.Ring.*)++psuedoScalar :: forall (p::Nat) (q::Nat) f. (Ord f, Algebra.Field.C f, SingI p, SingI q, SingI (p+q)) => Multivector p q f+psuedoScalar = one `e` [1..(toNatural d)] where+ d = fromIntegral (fromSing (sing :: Sing (p+q)) )::Word+ \end{code} Clifford numbers have a magnitude and absolute value:@@ -208,18 +213,18 @@ --magnitude :: (Algebra.Algebraic.C f) => Multivector f -> f magnitude = sqrt . compensatedSum' . map (\b -> (bScale b)^ 2) . mvTerms -instance (Algebra.Absolute.C f, Algebra.Algebraic.C f, Ord f, SingI n) => Algebra.Absolute.C (Multivector n f) where+instance (Algebra.Absolute.C f, Algebra.Algebraic.C f, Ord f, SingI p, SingI q) => Algebra.Absolute.C (Multivector p q f) where abs v = magnitude v `e` [] signum (BladeSum [Blade scale []]) = scalar $ signum scale signum (BladeSum []) = scalar Algebra.Additive.zero -instance (Algebra.Field.C f, Ord f, SingI n) => Algebra.Module.C f (Multivector n f) where+instance (Algebra.Field.C f, Ord f, SingI p, SingI q) => Algebra.Module.C f (Multivector p q f) where -- (*>) zero v = Algebra.Additive.zero (*>) s v = v & mvTerms & map (bladeScaleLeft s) & BladeSum -(/) :: (Algebra.Field.C f, Ord f, SingI n) => Multivector n f -> f -> Multivector n f+(/) :: (Algebra.Field.C f, Ord f, SingI p, SingI q) => Multivector p q f -> f -> Multivector p q f (/) v d = BladeSum $ map (bladeScaleLeft (NPN.recip d)) $ mvTerms v --Algebra.Field.recip d *> v (</) n d = Numeric.Clifford.Multivector.inverse d * n@@ -306,7 +311,7 @@ inverse a = assert (a /= zero) $ reverseMultivector a Numeric.Clifford.Multivector./ bScale (head $ mvTerms (a * reverseMultivector a)) recip=Numeric.Clifford.Multivector.inverse -instance (Algebra.Field.C f, Ord f, SingI n) => Algebra.OccasionallyScalar.C f (Multivector n f) where+instance (Algebra.Field.C f, Ord f, SingI p, SingI q) => Algebra.OccasionallyScalar.C f (Multivector p q f) where toScalar = bScale . bladeGetGrade 0 . head . mvTerms toMaybeScalar (BladeSum [Blade s []]) = Just s toMaybeScalar (BladeSum []) = Just Algebra.Additive.zero@@ -317,7 +322,7 @@ Also, we may as well implement the standard prelude Num interface. \begin{code}-instance (Algebra.Algebraic.C f, SingI n, Ord f) => PNum.Num (Multivector n f) where+instance (Algebra.Algebraic.C f, SingI p, SingI q, Ord f) => PNum.Num (Multivector p q f) where (+) = (Algebra.Additive.+) (-) = (Algebra.Additive.-) (*) = (Algebra.Ring.*)@@ -332,23 +337,23 @@ Let's use Newton or Halley iteration to find the principal n-th root :3 \begin{code}-root :: (Show f, Ord f, Algebra.Algebraic.C f, SingI d) => NPN.Integer -> Multivector d f -> Multivector d f+root :: (Show f, Ord f, Algebra.Algebraic.C f, SingI p, SingI q) => NPN.Integer -> Multivector p q f -> Multivector p q f root n (BladeSum [Blade s []]) = scalar $ Algebra.Algebraic.root n s root n a@(BladeSum _) = converge $ rootIterationsStart n a one -rootIterationsStart ::(Ord f, Show f, Algebra.Algebraic.C f)=> NPN.Integer -> Multivector d f -> Multivector d f -> [Multivector d f]+rootIterationsStart ::(Ord f, Show f, Algebra.Algebraic.C f)=> NPN.Integer -> Multivector p q f -> Multivector p q f -> [Multivector p q f] rootIterationsStart n a@(BladeSum (Blade s [] :xs)) one = rootHalleysIterations n a g where g = if s >= NPN.zero then one else Algebra.Ring.one `e` [1,2] --BladeSum[Blade Algebra.Ring.one [1,2]] rootIterationsStart n a@(BladeSum _) g = rootHalleysIterations n a g -rootNewtonIterations :: (Algebra.Field.C f, Ord f, SingI d) => NPN.Integer -> Multivector d f -> Multivector d f -> [Multivector d f]+rootNewtonIterations :: (Algebra.Field.C f, Ord f, SingI p, SingI q) => NPN.Integer -> Multivector p q f -> Multivector p q f -> [Multivector p q f] rootNewtonIterations n a = iterate xkplus1 where xkplus1 xk = xk + deltaxk xk deltaxk xk = oneOverN * ((Numeric.Clifford.Multivector.inverse (xk ^ (n - one))* a) - xk) oneOverN = scalar $ NPN.recip $ fromInteger n -rootHalleysIterations :: (Show a, Ord a, Algebra.Algebraic.C a, SingI d) => NPN.Integer -> Multivector d a -> Multivector d a -> [Multivector d a]+rootHalleysIterations :: (Show a, Ord a, Algebra.Algebraic.C a, SingI p, SingI q) => NPN.Integer -> Multivector p q a -> Multivector p q a -> [Multivector p q a] rootHalleysIterations n a = halleysMethod f f' f'' where f x = a - (x^n) f' x = fromInteger (-n) * (x^(n-1))@@ -361,7 +366,7 @@ down = denominator ratio -halleysMethod :: (Show a, Ord a, Algebra.Algebraic.C a, SingI d) => (Multivector d a -> Multivector d a) -> (Multivector d a -> Multivector d a) -> (Multivector d a -> Multivector d a) -> Multivector d a -> [Multivector d a]+halleysMethod :: (Show a, Ord a, Algebra.Algebraic.C a, SingI p, SingI q) => (Multivector p q a -> Multivector p q a) -> (Multivector p q a -> Multivector p q a) -> (Multivector p q a -> Multivector p q a) -> Multivector p q a -> [Multivector p q a] halleysMethod f f' f'' = iterate update where update x = x - (numerator x * Numeric.Clifford.Multivector.inverse (denominator x)) where numerator x = multiplyList [2, fx, dfx]
src/Numeric/Clifford/NumericIntegration.lhs view
@@ -97,8 +97,8 @@ showOutput name x = trace ("output of " ++ name ++" is " ++ show x) x -convergeTolLists :: (Ord f, Algebra.Absolute.C f, Algebra.Algebraic.C f, Show f, SingI d) - => f -> [[Multivector d f]] -> [Multivector d f]+convergeTolLists :: (Ord f, Algebra.Absolute.C f, Algebra.Algebraic.C f, Show f, SingI p, SingI q) + => f -> [[Multivector p q f]] -> [Multivector p q f] convergeTolLists t [] = error "converge: empty list" convergeTolLists t xs = fromMaybe empty (convergeBy check Just xs) where@@ -110,17 +110,17 @@ magnitude (sumList $ (zipWith (\x y -> NPN.abs (x-y)) b c))) <= t) = showOutput ("convergence with tolerance "++ show t )$ Just c check _ = Nothing -type RKStepper (d::Nat) t stateType = - (Ord t, Show t, Algebra.Module.C t (Multivector d t), Algebra.Field.C t) => +type RKStepper (p::Nat) (q::Nat) t stateType = + (Ord t, Show t, Algebra.Module.C t (Multivector p q t), Algebra.Field.C t) => t -> (t -> stateType -> stateType) -> - ([Multivector d t] -> stateType) -> - (stateType ->[Multivector d t]) ->+ ([Multivector p q t] -> stateType) -> + (stateType ->[Multivector p q t]) -> (t,stateType) -> (t,stateType) data ButcherTableau f = ButcherTableau {_tableauA :: [[f]], _tableauB :: [f], _tableauC :: [f]} makeLenses ''ButcherTableau -type ConvergerFunction f = forall (d::Nat) f . [[Multivector d f]] -> [Multivector d f]+type ConvergerFunction f = forall (p::Nat) (q::Nat) f . [[Multivector p q f]] -> [Multivector p q f] type AdaptiveStepSizeFunction f state = f -> state -> f data RKAttribute f state = Explicit@@ -135,9 +135,9 @@ $( derive makeIs ''RKAttribute) -genericRKMethod :: forall (d::Nat) t stateType . - ( Ord t, Show t, Algebra.Module.C t (Multivector d t),Algebra.Absolute.C t, Algebra.Algebraic.C t, SingI d)- => ButcherTableau t -> [RKAttribute t stateType] -> RKStepper d t stateType+genericRKMethod :: forall (p::Nat) (q::Nat) t stateType . + ( Ord t, Show t, Algebra.Module.C t (Multivector p q t),Algebra.Absolute.C t, Algebra.Algebraic.C t, SingI p, SingI q)+ => ButcherTableau t -> [RKAttribute t stateType] -> RKStepper p q t stateType genericRKMethod tableau attributes = rkMethodImplicitFixedPoint where s :: Int s = length (_tableauC tableau)@@ -151,10 +151,10 @@ b i = l !! (i - 1) where l = _tableauB tableau - sumListOfLists :: [[Multivector d t]] -> [Multivector d t]+ sumListOfLists :: [[Multivector p q t]] -> [Multivector p q t] sumListOfLists = map sumList . transpose - converger :: [[Multivector d t]] -> [Multivector d t]+ converger :: [[Multivector p q t]] -> [Multivector p q t] converger = case find (\x -> isConvergenceTolerance x || isConvergenceFunction x) attributes of Just (ConvergenceFunction conv) -> conv Just (ConvergenceTolerance tol) -> convergeTolLists (trace ("Convergence tolerance set to " ++ show tol)tol)@@ -165,13 +165,13 @@ Just (AdaptiveStepSize sigma) -> sigma Nothing -> (\_ _ -> one) - rkMethodImplicitFixedPoint :: RKStepper d t stateType+ rkMethodImplicitFixedPoint :: RKStepper p q t stateType rkMethodImplicitFixedPoint h f project unproject (time, state) = (time + (stepSizeAdapter time state)*h*(c s), newState) where- zi :: Int -> [Multivector d t]+ zi :: Int -> [Multivector p q t] zi i = (\out -> trace ("initialGuess is " ++ show initialGuess++" whereas the final one is " ++ show out) out) $ assert (i <= s && i>= 1) $ converger $ iterate (zkp1 i) initialGuess where- initialGuess :: [Multivector d t]+ initialGuess :: [Multivector p q t] initialGuess = if i == 1 || null (zi (i-1)) then map (h'*>) $ unproject $ f guessTime state else zi (i-1) adaptiveStepSizeFraction :: t adaptiveStepSizeFraction = stepSizeAdapter time state@@ -179,18 +179,18 @@ h' = adaptiveStepSizeFraction * h * (c i) guessTime :: t guessTime = time + h'- zkp1 :: NPN.Int -> [Multivector d t] -> [Multivector d t]+ zkp1 :: NPN.Int -> [Multivector p q t] -> [Multivector p q t] zkp1 i zk = map (h*>) (sumOfJs i zk) where- sumOfJs :: Int -> [Multivector d t] -> [Multivector d t]+ sumOfJs :: Int -> [Multivector p q t] -> [Multivector p q t] sumOfJs i zk = sumListOfLists $ map (scaledByAij zk) (a i) where - scaledByAij :: [Multivector d t] -> t -> [Multivector d t]+ scaledByAij :: [Multivector p q t] -> t -> [Multivector p q t] scaledByAij guess a = map (a*>) $ evalDerivatives guessTime $ elementAdd state' guess - state' = unproject state :: [Multivector d t]+ state' = unproject state :: [Multivector p q t] newState :: stateType newState = project $ elementAdd state' (assert (not $ null dy) dy)- dy = sumListOfLists [map ((b i) *>) (zi i) | i <- [1..s]] :: [Multivector d t]- evalDerivatives :: t -> [Multivector d t] -> [Multivector d t]+ dy = sumListOfLists [map ((b i) *>) (zi i) | i <- [1..s]] :: [Multivector p q t]+ evalDerivatives :: t -> [Multivector p q t] -> [Multivector p q t] evalDerivatives time stateAtTime= unproject $ (f time) $ project stateAtTime
src/Numeric/Clifford/NumericIntegration/DefaultIntegrators.hs view
@@ -36,6 +36,11 @@ import Test.QuickCheck --trace _ a = a +gaussLegendreFourthOrderTableau = ButcherTableau [[0.25::NPN.Double, 0.25 - (1.0 NPN./6)* sqrt 3], [0.25 + (1.0 NPN./ 6) * sqrt 3, 0.25]] [0.5, 0.5] [0.5 - (1.0 NPN./6)* sqrt 3, 0.5 + (1.0 NPN./ 6)* sqrt 3]+gaussLegendreFourthOrder h f (t, state) = impl h f id id (t,state) where+ impl= genericRKMethod gaussLegendreFourthOrderTableau []++ rk4ClassicalFromTableau h f (t,state) = impl h f id id (t, state) where impl = genericRKMethod rk4ClassicalTableau [] implicitEulerMethod h f (t, state) = impl h f id id (t, state) where@@ -55,7 +60,7 @@ lobattoIIIBFourthOrder h f (t, state) = impl h f id id (t, state) where impl = genericRKMethod lobattoIIIBFourthOrderTableau [] -rk4Classical :: (Ord a, Algebra.Algebraic.C a, SingI d) => stateType -> a -> (stateType->stateType) -> ([Multivector d a] -> stateType) -> (stateType -> [Multivector d a]) -> stateType+rk4Classical :: (Ord a, Algebra.Algebraic.C a, SingI p, SingI q) => stateType -> a -> (stateType->stateType) -> ([Multivector p q a] -> stateType) -> (stateType -> [Multivector p q a]) -> stateType rk4Classical state h f project unproject = project newState where update = map (\(k1', k2', k3', k4') -> sumList [k1',2*k2',2*k3',k4'] MV./ Algebra.Ring.fromInteger 6) $ zip4 k1 k2 k3 k4 newState = zipWith (+) state' update