diff --git a/learn-physics.cabal b/learn-physics.cabal --- a/learn-physics.cabal +++ b/learn-physics.cabal @@ -1,5 +1,5 @@ Name: learn-physics -Version: 0.6.0.0 +Version: 0.6.0.1 Synopsis: Haskell code for learning physics Description: A library of functions for vector calculus, calculation of electric field, electric flux, @@ -12,7 +12,7 @@ Category: Physics Build-type: Simple Cabal-version: >=1.8 -Tested-with: GHC == 7.10.2 +Tested-with: GHC >= 7.10.2 && <= 7.10.3 Library Exposed-modules: Physics.Learn.Charge Physics.Learn.Current @@ -39,7 +39,7 @@ Physics.Learn.Visual.PlotTools Physics.Learn.Visual.VisTools Physics.Learn.Visual.GlossTools - Build-depends: base >= 4.7 && < 4.9, + Build-depends: base >= 4.7 && < 4.10, vector-space >= 0.8.4 && < 0.11, not-gloss >= 0.5.0.4 && < 0.8, spatial-math >= 0.1.7 && < 0.3, @@ -56,44 +56,44 @@ Executable learn-physics-PlaneWave Main-is: examples/src/PlaneWave.hs Build-depends: not-gloss >= 0.7.4 && < 0.8, - base >= 4.5 && < 4.9, + base >= 4.5 && < 4.10, learn-physics Executable learn-physics-eFieldLine3D Main-is: examples/src/eFieldLine3D.hs Build-depends: not-gloss >= 0.7.4 && < 0.8, - base >= 4.5 && < 4.9, + base >= 4.5 && < 4.10, learn-physics Executable learn-physics-LorentzForceSimulation Main-is: examples/src/LorentzForceSimulation.hs Build-depends: not-gloss >= 0.7.4 && < 0.8, spatial-math >= 0.2 && < 0.3, - base >= 4.5 && < 4.9, + base >= 4.5 && < 4.10, learn-physics Executable learn-physics-BCircularLoop Main-is: examples/src/BCircularLoop.hs Build-depends: not-gloss >= 0.7.4 && < 0.8, - base >= 4.5 && < 4.9, + base >= 4.5 && < 4.10, learn-physics Executable learn-physics-sunEarth Main-is: examples/src/sunEarthRK4.hs Build-depends: gloss >= 1.8, - base >= 4.5 && < 4.9, + base >= 4.5 && < 4.10, learn-physics Executable learn-physics-eFieldLine2D Main-is: examples/src/eFieldLine2D.hs Build-depends: gloss >= 1.8, - base >= 4.5 && < 4.9, + base >= 4.5 && < 4.10, learn-physics Executable learn-physics-Projectile Main-is: examples/src/Projectile.hs Build-depends: gnuplot >= 0.5 && < 0.6, - base >= 4.5 && < 4.9, + base >= 4.5 && < 4.10, learn-physics Executable learn-physics-NMR diff --git a/src/Physics/Learn/BlochSphere.hs b/src/Physics/Learn/BlochSphere.hs --- a/src/Physics/Learn/BlochSphere.hs +++ b/src/Physics/Learn/BlochSphere.hs @@ -15,24 +15,27 @@ module Physics.Learn.BlochSphere ( VisObj , toPos + , ketToPos , staticBlochSphere , displayStaticState , animatedBlochSphere , simulateBlochSphere + , simulateBlochSphereK , stateProp + , statePropK , evolutionBlochSphere + , evolutionBlochSphereK , hamRabi ) where -import Physics.Learn.QuantumMat - ( sy - , sz - , xp - , yp - , ym - , zm - , timeEv +import qualified Physics.Learn.QuantumMat as M +import qualified Physics.Learn.Ket as K +import Physics.Learn.Ket + ( Ket + , Operator + , (<>) + , dagger ) import Numeric.LinearAlgebra ( Vector @@ -40,7 +43,7 @@ , C , iC -- , (<>) -- matrix multiplication - , (|>) -- vector definition +-- , (|>) -- vector definition , (!) -- vector element access , (><) -- matrix definition , scale @@ -72,6 +75,22 @@ , red ) +{- +3 ways to specify the state of a spin-1/2 particle: +Vector C +Ket +Position (Bloch vector) + +2 ways to specify a Hamiltonian: +Matrix C +Operator + +3 choices for Vis' world: +(Float, Vector C) +(Float, Ket) +(Float, Position) +-} + -- | A Vis object. type VisObj = VisObject Double @@ -87,6 +106,18 @@ (2 * imagPart (conjugate z1 * z2)) (realPart (conjugate z1 * z1 - conjugate z2 * z2)) +-- | Convert a qubit ket +-- into Bloch (x,y,z) coordinates. +ketToPos :: Ket -> Position +ketToPos psi + = if K.dim psi /= 2 + then error "ketToPos only for qubit kets" + else let z1 = dagger K.zp <> psi + z2 = dagger K.zm <> psi + in cart (2 * realPart (conjugate z1 * z2)) + (2 * imagPart (conjugate z1 * z2)) + (realPart (conjugate z1 * z1 - conjugate z2 * z2)) + -- | A static 'VisObj' for the state of a qubit. staticBlochSphere :: Position -> VisObj staticBlochSphere r @@ -103,15 +134,6 @@ displayStaticState :: Vector C -> IO () displayStaticState = displayStaticBlochSphere . toPos -displayxp :: IO () -displayxp = displayStaticState xp - -displayyp :: IO () -displayyp = displayStaticState yp - -displayym :: IO () -displayym = displayStaticState ym - -- | Given a Bloch vector as a function of time, -- return a 'VisObj' as a function of time. animatedBlochSphere :: (Double -> Position) -> (Float -> VisObj) @@ -126,6 +148,14 @@ simulateBlochSphere sampleRate initial statePropFunc = simulate myOptions sampleRate (0,initial) (staticBlochSphere . toPos . snd) statePropFunc +-- | Given a sample rate, initial qubit state ket, and +-- state propagation function, produce a simulation. +-- The 'Float' in the state propagation function is the time +-- since the beginning of the simulation. +simulateBlochSphereK :: Double -> Ket -> (Float -> (Float,Ket) -> (Float,Ket)) -> IO () +simulateBlochSphereK sampleRate initial statePropFuncK + = simulate myOptions sampleRate (0,initial) (staticBlochSphere . ketToPos . snd) statePropFuncK + {- -- | Given a sample rate, initial qubit state vector, and -- state propagation function, produce a simulation. @@ -139,29 +169,45 @@ -- | Produce a state propagation function from a time-dependent Hamiltonian. stateProp :: (Double -> Matrix C) -> Float -> (Float,Vector C) -> (Float,Vector C) stateProp ham tNew (tOld,v) - = (tNew, timeEv (realToFrac dt) (ham tMid) v) + = (tNew, M.timeEv (realToFrac dt) (ham tMid) v) where dt = tNew - tOld tMid = realToFrac $ (tNew + tOld) / 2 +-- | Produce a state propagation function from a time-dependent Hamiltonian. +statePropK :: (Double -> Operator) -> Float -> (Float,Ket) -> (Float,Ket) +statePropK ham tNew (tOld,psi) + = (tNew, K.timeEv (realToFrac dt) (ham tMid) psi) + where + dt = tNew - tOld + tMid = realToFrac $ (tNew + tOld) / 2 + -- | Given an initial qubit state and a time-dependent Hamiltonian, -- produce a visualization. evolutionBlochSphere :: Vector C -> (Double -> Matrix C) -> IO () evolutionBlochSphere psi0 ham = simulateBlochSphere 0.01 psi0 (stateProp ham) +-- | Given an initial qubit ket and a time-dependent Hamiltonian, +-- produce a visualization. +evolutionBlochSphereK :: Ket -> (Double -> Operator) -> IO () +evolutionBlochSphereK psi0 ham + = simulateBlochSphereK 0.01 psi0 (statePropK ham) + myOptions :: Options myOptions = defaultOpts {optWindowName = "Bloch Sphere" ,optInitialCamera = Just (Camera0 75 20 4)} +{- staticBz1 :: IO () -staticBz1 = evolutionBlochSphere xp (const (scale 0.9 sz)) +staticBz1 = evolutionBlochSphere M.xp (const (scale 0.9 M.sz)) staticBz2 :: IO () -staticBz2 = evolutionBlochSphere ((2|>) [(cos (pi / 8)), (sin (pi / 8))]) (const sz) +staticBz2 = evolutionBlochSphere ((2|>) [(cos (pi / 8)), (sin (pi / 8))]) (const M.sz) staticBy1 :: IO () -staticBy1 = evolutionBlochSphere xp (const sy) +staticBy1 = evolutionBlochSphere M.xp (const M.sy) +-} -- | Hamiltonian for nuclear magnetic resonance. -- Explain omega0, omegaR, omega. diff --git a/src/Physics/Learn/Ket.hs b/src/Physics/Learn/Ket.hs --- a/src/Physics/Learn/Ket.hs +++ b/src/Physics/Learn/Ket.hs @@ -17,16 +17,15 @@ -- a Ket layer on top of QuantumMat module Physics.Learn.Ket - ( Ket + ( + -- * Basic data types + C + , i + , magnitude + , Ket , Bra , Operator - , Mult(..) - , Dagger(..) - , Representable(..) - , OrthonormalBasis - , makeOB - , listBasis - , size + -- * Kets for spin-1/2 particles , xp , xm , yp @@ -35,14 +34,41 @@ , zm , np , nm - , xBasis - , yBasis - , zBasis + -- * Operators for spin-1/2 particles , sx , sy , sz - , prob - , probs + , sn + , sn' + -- * Quantum Dynamics + , timeEvOp + , timeEv + -- * Composition + , Kron(..) + -- * Measurement + , possibleOutcomes + , outcomesProjectors + , outcomesProbabilities +-- , prob +-- , probs + -- * Generic multiplication + , Mult(..) + -- * Adjoint operation + , Dagger(..) + -- * Normalization + , HasNorm(..) + -- * Representation + , Representable(..) + -- * Orthonormal bases + , OrthonormalBasis + , makeOB + , listBasis + , size + -- * Orthonormal bases for spin-1/2 particles + , xBasis + , yBasis + , zBasis + , nBasis -- , angularMomentumXMatrix -- , angularMomentumYMatrix -- , angularMomentumZMatrix @@ -60,9 +86,9 @@ -- We try to import only from QuantumMat -- and not from Numeric.LinearAlgebra +import qualified Data.Complex as C import Data.Complex ( Complex(..) - , magnitude , conjugate ) import qualified Physics.Learn.QuantumMat as M @@ -72,7 +98,6 @@ , Matrix , (#>) , (<#) - , couter , conjugateTranspose , scaleV , scaleM @@ -100,7 +125,12 @@ data Operator = Operator (Matrix C) instance Show Operator where - show _ = "\nTry 'rep zBasis '" + show op = + let message = "Use 'rep '." + in if dim op == 2 + then "Representation in zBasis:\n" ++ + show (rep zBasis op) ++ "\n" ++ message + else message -- | A bra vector describes the state of a quantum system. data Bra = Bra (Vector C) @@ -108,6 +138,12 @@ instance Show Bra where show _ = "\nTry 'rep zBasis '" +magnitude :: C -> Double +magnitude = C.magnitude + +i :: C +i = 0 :+ 1 + -- | Generic multiplication including inner product, -- outer product, operator product, and whatever else makes sense. -- No conjugation takes place in this operation. @@ -148,7 +184,9 @@ = Ket (matrixOp #> matrixKet) instance Mult Ket Bra Operator where - Ket k <> Bra b = Operator (couter k b) + Ket k <> Bra b = + Operator + (fromLists [[ x*y | y <- toList b] | x <- toList k]) instance Mult Operator Operator Operator where Operator m1 <> Operator m2 = Operator (m1 M.<> m2) @@ -209,17 +247,6 @@ norm (Bra v) = M.norm v normalize b = (1 / norm b :+ 0) <> b -{- -class HasDim a where - dim :: a -> Int - -instance HasDim Ket where - dim (Ket v) = M.dim v - -instance HasDim Bra where - dim (Bra v) = M.dim v --} - -- | An orthonormal basis of kets. newtype OrthonormalBasis = OB [Ket] deriving (Show) @@ -256,14 +283,6 @@ dim (Operator m) = let (p,q) = M.size m in if p == q then p else error "dim: non-square operator" -prob :: Ket -> Ket -> Double -prob k1 k2 = magnitude c ** 2 - where - c = dagger k1 <> k2 - -probs :: OrthonormalBasis -> Ket -> [Double] -probs (OB ks) k = map (\bk -> let c = dagger bk <> k in magnitude c ** 2) ks - -------------- -- Spin 1/2 -- -------------- @@ -314,15 +333,23 @@ = (sin (theta / 2) :+ 0) <> zp - (cos (theta / 2) :+ 0) * (cos phi :+ sin phi) <> zm +-- | The orthonormal basis composed of 'xp' and 'xm'. xBasis :: OrthonormalBasis xBasis = makeOB [xp,xm] +-- | The orthonormal basis composed of 'yp' and 'ym'. yBasis :: OrthonormalBasis yBasis = makeOB [yp,ym] +-- | The orthonormal basis composed of 'zp' and 'zm'. zBasis :: OrthonormalBasis zBasis = makeOB [zp,zm] +-- | Given spherical polar angle theta and azimuthal angle phi, +-- the orthonormal basis composed of 'np' theta phi and 'nm' theta phi. +nBasis :: Double -> Double -> OrthonormalBasis +nBasis theta phi = makeOB [np theta phi,nm theta phi] + -- | The Pauli X operator. sx :: Operator sx = xp <> dagger xp - xm <> dagger xm @@ -334,6 +361,93 @@ -- | The Pauli Z operator. sz :: Operator sz = zp <> dagger zp - zm <> dagger zm + +-- | Pauli operator for an arbitrary direction given +-- by spherical coordinates theta and phi. +sn :: Double -> Double -> Operator +sn theta phi + = (sin theta * cos phi :+ 0) <> sx + + (sin theta * sin phi :+ 0) <> sy + + (cos theta :+ 0) <> sz + +-- | Alternative definition +-- of Pauli operator for an arbitrary direction. +sn' :: Double -> Double -> Operator +sn' theta phi + = np theta phi <> dagger (np theta phi) - + nm theta phi <> dagger (nm theta phi) + +---------------------- +-- Quantum Dynamics -- +---------------------- + +-- | Given a time step and a Hamiltonian operator, +-- produce a unitary time evolution operator. +-- Unless you really need the time evolution operator, +-- it is better to use 'timeEv', which gives the +-- same numerical results without doing an explicit +-- matrix inversion. The function assumes hbar = 1. +timeEvOp :: Double -> Operator -> Operator +timeEvOp dt (Operator m) = Operator (M.timeEvMat dt m) + +-- | Given a time step and a Hamiltonian operator, +-- advance the state ket using the Schrodinger equation. +-- This method should be faster than using 'timeEvOp' +-- since it solves a linear system rather than calculating +-- an inverse matrix. The function assumes hbar = 1. +timeEv :: Double -> Operator -> Ket -> Ket +timeEv dt (Operator m) (Ket k) = Ket $ M.timeEv dt m k + +----------------- +-- Composition -- +----------------- + +class Kron a where + kron :: a -> a -> a + +instance Kron Ket where + kron (Ket v1) (Ket v2) = Ket (M.kron v1 v2) + +instance Kron Bra where + kron (Bra v1) (Bra v2) = Bra (M.kron v1 v2) + +instance Kron Operator where + kron (Operator m1) (Operator m2) = Operator (M.kron m1 m2) + +----------------- +-- Measurement -- +----------------- + +-- | The possible outcomes of a measurement +-- of an observable. +-- These are the eigenvalues of the operator +-- of the observable. +possibleOutcomes :: Operator -> [Double] +possibleOutcomes (Operator observable) = M.possibleOutcomes observable + +-- | Given an obervable, return a list of pairs +-- of possible outcomes and projectors +-- for each outcome. +outcomesProjectors :: Operator -> [(Double,Operator)] +outcomesProjectors (Operator m) + = [(val,Operator p) | (val,p) <- M.outcomesProjectors m] + +-- | Given an observable and a state ket, return a list of pairs +-- of possible outcomes and probabilites +-- for each outcome. +outcomesProbabilities :: Operator -> Ket -> [(Double,Double)] +outcomesProbabilities (Operator m) (Ket v) + = M.outcomesProbabilities m v + +{- +prob :: Ket -> Ket -> Double +prob k1 k2 = magnitude c ** 2 + where + c = dagger k1 <> k2 + +probs :: OrthonormalBasis -> Ket -> [Double] +probs (OB ks) k = map (\bk -> let c = dagger bk <> k in magnitude c ** 2) ks +-} {- ---------------------------------------- diff --git a/src/Physics/Learn/QuantumMat.hs b/src/Physics/Learn/QuantumMat.hs --- a/src/Physics/Learn/QuantumMat.hs +++ b/src/Physics/Learn/QuantumMat.hs @@ -49,6 +49,7 @@ , fromLists , toLists , size + , matrixFunction -- * Density matrices , couter , dm @@ -56,10 +57,15 @@ , normalizeDM , oneQubitMixed -- * Quantum Dynamics - , timeEv , timeEvMat + , timeEv + , timeEvMatSpec + -- * Composition + , Kronecker(..) -- * Measurement , possibleOutcomes + , outcomesProjectors + , outcomesProbabilities -- * Vector and Matrix , Vector , Matrix @@ -70,6 +76,7 @@ ( C , Vector , Matrix + , Herm , iC -- square root of negative one , (><) -- matrix definition , ident @@ -79,6 +86,7 @@ , (<\>) , sym , eigenvaluesSH + , eigSH , cmap , takeDiag , conj @@ -94,6 +102,7 @@ import Data.Complex ( Complex(..) , magnitude + , realPart ) -- | The state resulting from a measurement of @@ -195,6 +204,10 @@ toList :: Vector C -> [C] toList = H.toList +-------------- +-- Matrices -- +-------------- + -- | The Pauli X matrix. sx :: Matrix C sx = (2><2) [ 0, 1 @@ -242,6 +255,15 @@ size :: Matrix C -> (Int,Int) size = H.size +-- | Apply a function to a matrix. +-- Assumes the matrix is a normal matrix (a matrix +-- with an orthonormal basis of eigenvectors). +matrixFunction :: (C -> C) -> Matrix C -> Matrix C +matrixFunction f m + = let (valv,vecm) = H.eig m + fvalv = fromList [f val | val <- toList valv] + in vecm <> H.diag fvalv <> tr vecm + ---------------------- -- Density Matrices -- ---------------------- @@ -274,7 +296,7 @@ -- produce a unitary time evolution matrix. -- Unless you really need the time evolution matrix, -- it is better to use 'timeEv', which gives the --- same numerical results with doing an explicit +-- same numerical results without doing an explicit -- matrix inversion. The function assumes hbar = 1. timeEvMat :: Double -> Matrix C -> Matrix C timeEvMat dt h @@ -297,7 +319,26 @@ identity = ident n in (identity + ah) <\> ((identity - ah) #> v) +-- | Given a Hamiltonian matrix, return a function from time +-- to evolution matrix. Uses spectral decomposition. +-- Assumes hbar = 1. +timeEvMatSpec :: Matrix C -> Double -> Matrix C +timeEvMatSpec m t = matrixFunction (\h -> exp(-iC * h * (t :+ 0))) m + ----------------- +-- Composition -- +----------------- + +class Kronecker a where + kron :: a -> a -> a + +instance H.Product t => Kronecker (Vector t) where + kron v1 v2 = H.fromList [c1 * c2 | c1 <- H.toList v1, c2 <- H.toList v2] + +instance H.Product t => Kronecker (Matrix t) where + kron = H.kronecker + +----------------- -- Measurement -- ----------------- @@ -308,6 +349,38 @@ possibleOutcomes :: Matrix C -> [Double] possibleOutcomes observable = H.toList $ eigenvaluesSH (sym observable) + +-- From a Hermitian matrix, a list of pairs of eigenvalues and eigenvectors. +valsVecs :: Herm C -> [(Double,Vector C)] +valsVecs h = let (valv,m) = eigSH h + vals = H.toList valv + vecs = map (conjV . fromList) $ toLists (conjugateTranspose m) + in zip vals vecs + +-- From a Hermitian matrix, a list of pairs of eigenvalues and projectors. +valsPs :: Herm C -> [(Double,Matrix C)] +valsPs h = [(val,couter vec vec) | (val,vec) <- valsVecs h] + +combineFst :: (Eq a, Num b) => [(a,b)] -> [(a,b)] +combineFst [] = [] +combineFst [p] = [p] +combineFst ((x1,m1):(x2,m2):ps) + = if x1 == x2 + then combineFst ((x1,m1+m2):ps) + else (x1,m1):combineFst ((x2,m2):ps) + +-- | Given an obervable, return a list of pairs +-- of possible outcomes and projectors +-- for each outcome. +outcomesProjectors :: Matrix C -> [(Double,Matrix C)] +outcomesProjectors m = combineFst (valsPs (sym m)) + +-- | Given an observable and a state vector, return a list of pairs +-- of possible outcomes and probabilites +-- for each outcome. +outcomesProbabilities :: Matrix C -> Vector C -> [(Double,Double)] +outcomesProbabilities m v + = [(a,realPart (inner v (p #> v))) | (a,p) <- outcomesProjectors m] ------------------ -- Gram-Schmidt --