From aef0333b5180ea79e539bd53194f1dfed20b7db5 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Sun, 7 Feb 2010 09:29:50 +0000 Subject: added odeSolve --- examples/HMatrixReal.hs | 104 ------------------------------------------------ examples/Real.hs | 104 ++++++++++++++++++++++++++++++++++++++++++++++++ examples/ode.hs | 34 ++++++++++++++++ 3 files changed, 138 insertions(+), 104 deletions(-) delete mode 100644 examples/HMatrixReal.hs create mode 100644 examples/Real.hs create mode 100644 examples/ode.hs (limited to 'examples') diff --git a/examples/HMatrixReal.hs b/examples/HMatrixReal.hs deleted file mode 100644 index 0edff10..0000000 --- a/examples/HMatrixReal.hs +++ /dev/null @@ -1,104 +0,0 @@ - --- Alternative interface and utilities for creation of real arrays, useful to work in interactive mode. - -module HMatrixReal( - module Numeric.LinearAlgebra, - (<>), (*>), (<*), (<\>), (\>), - vector, - eye, - zeros, ones, - diagl, - row, - col, - (#),(&), (//), blocks, - rand -) where - -import Numeric.LinearAlgebra hiding ((<>), (<|>), (<->), (<\>), (.*), (*/)) -import System.Random(randomIO) - -infixl 7 <> --- | Matrix product ('multiply') -(<>) :: Field t => Matrix t -> Matrix t -> Matrix t -(<>) = multiply - -infixl 7 *> --- | matrix x vector -(*>) :: Field t => Matrix t -> Vector t -> Vector t -m *> v = flatten $ m <> (asColumn v) - -infixl 7 <* --- | vector x matrix -(<*) :: Field t => Vector t -> Matrix t -> Vector t -v <* m = flatten $ (asRow v) <> m - - --- | Least squares solution of a linear system for several right-hand sides, similar to the \\ operator of Matlab\/Octave. (\<\\\>) = 'linearSolveSVD'. -(<\>) :: (Field a) => Matrix a -> Matrix a -> Matrix a -infixl 7 <\> -(<\>) = linearSolveSVD - --- | Least squares solution of a linear system for a single right-hand side. See '(\<\\\>)'. -(\>) :: (Field a) => Matrix a -> Vector a -> Vector a -infixl 7 \> -m \> v = flatten (m <\> reshape 1 v) - --- | Pseudorandom matrix with uniform elements between 0 and 1. -rand :: Int -- ^ rows - -> Int -- ^ columns - -> IO (Matrix Double) -rand r c = do - seed <- randomIO - return (reshape c $ randomVector seed Uniform (r*c)) - --- | Real identity matrix. -eye :: Int -> Matrix Double -eye = ident - --- | Create a real vector from a list. -vector :: [Double] -> Vector Double -vector = fromList - --- | Create a real diagonal matrix from a list. -diagl :: [Double] -> Matrix Double -diagl = diag . vector - --- | Create a matrix or zeros. -zeros :: Int -- ^ rows - -> Int -- ^ columns - -> Matrix Double -zeros r c = reshape c (constant 0 (r*c)) - --- | Create a matrix or ones. -ones :: Int -- ^ rows - -> Int -- ^ columns - -> Matrix Double -ones r c = reshape c (constant 1 (r*c)) - --- | Concatenation of real vectors. -infixl 9 # -(#) :: Vector Double -> Vector Double -> Vector Double -a # b = join [a,b] - --- | Horizontal concatenation of real matrices. -infixl 8 & -(&) :: Matrix Double -> Matrix Double -> Matrix Double -a & b = fromBlocks [[a,b]] - --- | Vertical concatenation of real matrices. -infixl 7 // -(//) :: Matrix Double -> Matrix Double -> Matrix Double -a // b = fromBlocks [[a],[b]] - --- | Real block matrix from a rectangular list of lists. -blocks :: [[Matrix Double]] -> Matrix Double -blocks = fromBlocks - --- | A real matrix with a single row, create from a list of elements. -row :: [Double] -> Matrix Double -row = asRow . vector - --- | A real matrix with a single column, created from a list of elements. -col :: [Double] -> Matrix Double -col = asColumn . vector - diff --git a/examples/Real.hs b/examples/Real.hs new file mode 100644 index 0000000..b32a961 --- /dev/null +++ b/examples/Real.hs @@ -0,0 +1,104 @@ + +-- Alternative interface and utilities for creation of real arrays, useful to work in interactive mode. + +module Real( + module Numeric.LinearAlgebra, + (<>), (*>), (<*), (<\>), (\>), + vector, + eye, + zeros, ones, + diagl, + row, + col, + (#),(&), (//), blocks, + rand +) where + +import Numeric.LinearAlgebra hiding ((<>), (<|>), (<->), (<\>), (.*), (*/)) +import System.Random(randomIO) + +infixl 7 <> +-- | Matrix product ('multiply') +(<>) :: Field t => Matrix t -> Matrix t -> Matrix t +(<>) = multiply + +infixl 7 *> +-- | matrix x vector +(*>) :: Field t => Matrix t -> Vector t -> Vector t +m *> v = flatten $ m <> (asColumn v) + +infixl 7 <* +-- | vector x matrix +(<*) :: Field t => Vector t -> Matrix t -> Vector t +v <* m = flatten $ (asRow v) <> m + + +-- | Least squares solution of a linear system for several right-hand sides, similar to the \\ operator of Matlab\/Octave. (\<\\\>) = 'linearSolveSVD'. +(<\>) :: (Field a) => Matrix a -> Matrix a -> Matrix a +infixl 7 <\> +(<\>) = linearSolveSVD + +-- | Least squares solution of a linear system for a single right-hand side. See '(\<\\\>)'. +(\>) :: (Field a) => Matrix a -> Vector a -> Vector a +infixl 7 \> +m \> v = flatten (m <\> reshape 1 v) + +-- | Pseudorandom matrix with uniform elements between 0 and 1. +rand :: Int -- ^ rows + -> Int -- ^ columns + -> IO (Matrix Double) +rand r c = do + seed <- randomIO + return (reshape c $ randomVector seed Uniform (r*c)) + +-- | Real identity matrix. +eye :: Int -> Matrix Double +eye = ident + +-- | Create a real vector from a list. +vector :: [Double] -> Vector Double +vector = fromList + +-- | Create a real diagonal matrix from a list. +diagl :: [Double] -> Matrix Double +diagl = diag . vector + +-- | Create a matrix or zeros. +zeros :: Int -- ^ rows + -> Int -- ^ columns + -> Matrix Double +zeros r c = reshape c (constant 0 (r*c)) + +-- | Create a matrix or ones. +ones :: Int -- ^ rows + -> Int -- ^ columns + -> Matrix Double +ones r c = reshape c (constant 1 (r*c)) + +-- | Concatenation of real vectors. +infixl 9 # +(#) :: Vector Double -> Vector Double -> Vector Double +a # b = join [a,b] + +-- | Horizontal concatenation of real matrices. +infixl 8 & +(&) :: Matrix Double -> Matrix Double -> Matrix Double +a & b = fromBlocks [[a,b]] + +-- | Vertical concatenation of real matrices. +infixl 7 // +(//) :: Matrix Double -> Matrix Double -> Matrix Double +a // b = fromBlocks [[a],[b]] + +-- | Real block matrix from a rectangular list of lists. +blocks :: [[Matrix Double]] -> Matrix Double +blocks = fromBlocks + +-- | A real matrix with a single row, create from a list of elements. +row :: [Double] -> Matrix Double +row = asRow . vector + +-- | A real matrix with a single column, created from a list of elements. +col :: [Double] -> Matrix Double +col = asColumn . vector + diff --git a/examples/ode.hs b/examples/ode.hs new file mode 100644 index 0000000..082c46c --- /dev/null +++ b/examples/ode.hs @@ -0,0 +1,34 @@ +import Numeric.GSL.ODE +import Numeric.LinearAlgebra +import Graphics.Plot + +vanderpol mu = do + let xdot mu t [x,v] = [v, -x + mu * v * (1-x^2)] + ts = linspace 1000 (0,50) + sol = toColumns $ odeSolve (xdot mu) [1,0] ts + mplot (ts : sol) + mplot sol + + +harmonic w d = do + let xdot w d t [x,v] = [v, a*x + b*v] where a = -w^2; b = -2*d*w + ts = linspace 100 (0,20) + sol = odeSolve (xdot w d) [1,0] ts + mplot (ts : toColumns sol) + + +kepler v a = mplot (take 2 $ toColumns sol) where + xdot t [x,y,vx,vy] = [vx,vy,x*k,y*k] + where g=1 + k=(-g)*(x*x+y*y)**(-1.5) + ts = linspace 100 (0,30) + sol = odeSolve xdot [4, 0, v * cos (a*degree), v * sin (a*degree)] ts + degree = pi/180 + + +main = do + vanderpol 2 + harmonic 1 0 + harmonic 1 0.1 + kepler 0.3 60 + kepler 0.4 70 -- cgit v1.2.3