aurellem

Building a Classical Mechanics Library in Clojure

Table of Contents

1 Useful data types

1.1 Complex numbers

1.2 Power series

1.3 Tuples and tensors

1.3.1 Tuples are “sequences with spin”

(in-ns 'sicm.utils)

;; Let some objects have spin

(defprotocol Spinning
  (up? [this])
  (down? [this]))

(defn spin
  "Returns the spin of the Spinning s, either :up or :down"
  [#^Spinning s]
  (cond (up? s) :up (down? s) :down))


;; DEFINITION: A tuple is a sequence with spin

(deftype Tuple
  [spin coll]

  clojure.lang.Seqable
  (seq [this] (seq (.coll this)))

  clojure.lang.Counted
  (count [this] (count (.coll this)))

  Spinning
  (up? [this] (= ::up (.spin this)))
  (down? [this] (= ::down (.spin this))))

(defmethod print-method Tuple
  [o w]
  (print-simple
   (if (up? o)
     (str "u" (.coll o))
     (str "d" (vec(.coll o))))
   w))


;; CONSTRUCTORS

(defn up
  "Create a new up-tuple containing the contents of coll."
  [coll]
  (Tuple. ::up coll))       

(defn down
  "Create a new down-tuple containing the contents of coll."
  [coll]
  (Tuple. ::down coll))

1.3.2 Matrices

(in-ns 'sicm.utils)
(require 'incanter.core) ;; use incanter's fast matrices


(defn all-equal? [coll]
  (if (empty? (rest coll)) true
      (and (= (first coll) (second coll))
           (recur (rest coll)))))


(defprotocol Matrix
  (rows [matrix])
  (cols [matrix])
  (diagonal [matrix])
  (trace [matrix])
  (determinant [matrix])
  (transpose [matrix])
  (conjugate [matrix])
)

(extend-protocol Matrix incanter.Matrix
  (rows [rs] (map down (apply map vector (apply map vector rs))))
  (cols [rs] (map up (apply map vector rs)))
  (diagonal [matrix] (incanter.core/diag matrix) )
  (determinant [matrix] (incanter.core/det matrix))
  (trace [matrix] (incanter.core/trace matrix))
  (transpose [matrix] (incanter.core/trans matrix)))

(defn count-rows [matrix]
  ((comp count rows) matrix))

(defn count-cols [matrix]
  ((comp count cols) matrix))

(defn square? [matrix]
  (= (count-rows matrix) (count-cols matrix)))

(defn identity-matrix
  "Define a square matrix of size n-by-n with 1s along the diagonal and
  0s everywhere else."
  [n]
  (incanter.core/identity-matrix n))


(defn matrix-by-rows
  "Define a matrix by giving its rows."
  [& rows]
  (if
   (not (all-equal? (map count rows)))
   (throw (Exception. "All rows in a matrix must have the same number of elements."))
   (incanter.core/matrix (vec rows))))

(defn matrix-by-cols
  "Define a matrix by giving its columns"
  [& cols]
  (if (not (all-equal? (map count cols)))
   (throw (Exception. "All columns in a matrix must have the same number of elements."))
   (incanter.core/matrix (vec (apply map vector cols)))))

(defn identity-matrix
  "Define a square matrix of size n-by-n with 1s along the diagonal and
  0s everywhere else."
  [n]
  (incanter.core/identity-matrix n))

1.4 Generic Operations

(in-ns 'sicm.utils)
(use 'clojure.contrib.generic.arithmetic
     'clojure.contrib.generic.collection
     'clojure.contrib.generic.functor
     'clojure.contrib.generic.math-functions)

(defn numbers?
  "Returns true if all arguments are numbers, else false."
  [& xs]
  (every? number? xs))

(defn tuple-surgery
  "Applies the function f to the items of tuple and the additional
  arguments, if any. Returns a Tuple of the same type as tuple."
  [tuple f & xs]
  ((if (up? tuple) up down)
   (apply f (seq tuple) xs)))



;;; CONTRACTION collapses two compatible tuples into a number.

(defn contractible?
  "Returns true if the tuples a and b are compatible for contraction,
  else false. Tuples are compatible if they have the same number of
  components, they have opposite spins, and their elements are
  pairwise-compatible."
  [a b]
  (and
   (isa? (type a) Tuple)
   (isa? (type b) Tuple)
   (= (count a) (count b))
   (not= (spin a) (spin b))
   
   (not-any? false?
             (map #(or
                    (numbers? %1 %2)
                    (contractible? %1 %2))
                  a b))))

(defn contract
  "Contracts two tuples, returning the sum of the
  products of the corresponding items. Contraction is recursive on
  nested tuples."
  [a b]
  (if (not (contractible? a b))
    (throw
     (Exception. "Not compatible for contraction."))
    (reduce +
            (map
             (fn [x y]
               (if (numbers? x y)
                 (* x y)
                 (contract x y)))
             a b))))





(defmethod conj Tuple
  [tuple & xs]
  (tuple-surgery tuple #(apply conj % xs)))

(defmethod fmap Tuple
  [f tuple]
  (tuple-surgery tuple (partial map f)))



;; TODO: define Scalar, and add it to the hierarchy above Number and Complex

                                        
(defmethod * [Tuple Tuple]             ; tuple*tuple
  [a b]
  (if (contractible? a b)
    (contract a b)
    (map (partial * a) b)))


(defmethod * [java.lang.Number Tuple]  ;; scalar *  tuple
  [a x] (fmap (partial * a) x))

(defmethod * [Tuple java.lang.Number]
  [x a] (* a x))

(defmethod * [java.lang.Number incanter.Matrix] ;; scalar *  matrix
  [x M] (incanter.core/mult x M))

(defmethod * [incanter.Matrix java.lang.Number]
  [M x] (* x M))

(defmethod * [incanter.Matrix incanter.Matrix] ;; matrix * matrix
  [M1 M2]
  (incanter.core/mmult M1 M2))

(defmethod * [incanter.Matrix Tuple] ;; matrix * tuple
  [M v]
  (if (and (apply numbers? v) (up? v)) 
    (* M (matrix-by-cols v))
    (throw (Exception. "Currently, you can only multiply a matrix by a tuple of *numbers*"))
    ))

(defmethod * [Tuple incanter.Matrix] ;; tuple * Matrix
  [v M]
  (if (and (apply numbers? v) (down? v))
    (* (matrix-by-rows v) M)
    (throw (Exception. "Currently, you can only multiply a matrix by a tuple of *numbers*"))
    ))


(defmethod exp incanter.Matrix
  [M]
  (incanter.core/exp M))

2 Operators and Differentiation

2.1 Operators

2.2 Differential Terms and Sequences

(in-ns 'sicm.utils)
(use 'clojure.contrib.seq
     'clojure.contrib.generic.arithmetic
     'clojure.contrib.generic.collection
     'clojure.contrib.generic.math-functions)

;;

;; DEFINITION : Differential Term

;; A quantity with infinitesimal components, e.g. x, dxdy, 4dydz. The
;; coefficient of the quantity is returned by the 'coefficient' method,
;; while the sequence of differential parameters is returned by the
;; method 'partials'.

;; Instead of using (potentially ambiguous) letters to denote
;; differential parameters (dx,dy,dz), we use integers. So, dxdz becomes [0 2].

;; The coefficient can be any arithmetic object; the
;; partials must be a nonrepeating sorted sequence of nonnegative
;; integers.

(deftype DifferentialTerm [coefficient partials])

(defn differential-term
  "Make a differential term given pairs of coefficients and partials."
  [coefficient partials]
  (if (and (coll? partials) (every? #(and (integer? %) (not(neg? %))) partials)) 
    (DifferentialTerm. coefficient (set partials))
    (throw (java.lang.IllegalArgumentException. "Partials must be a collection of integers."))))


;; DEFINITION : Differential Sequence
;; A differential sequence is a sequence of differential terms, all with different partials.
;; Internally, it is a map from the partials of each term to their coefficients.

(deftype DifferentialSeq
  [terms]
  ;;clojure.lang.IPersistentMap
  clojure.lang.Associative
  (assoc [this key val]
    (DifferentialSeq.
     (cons (differential-term val key) terms)))
  (cons [this x]
        (DifferentialSeq. (cons x terms)))
  (containsKey [this key]
               (not(nil? (find-first #(= (.partials %) key) terms))))
  (count [this] (count (.terms this)))
  (empty [this] (DifferentialSeq. []))
  (entryAt [this key]
           ((juxt #(.partials %) #(.coefficient %))
            (find-first #(= (.partials %) key) terms)))
  (seq [this] (seq (.terms this))))


(defmethod fmap DifferentialSeq
  [f dseq]
  (DifferentialSeq. (fmap f (.terms dseq))))


;; BUILDING DIFFERENTIAL OBJECTS

(defn differential-seq
    "Define a differential sequence by specifying coefficient/partials
pairs, which are used to create differential terms to populate the
sequence."
  ([coefficient partials]
     (DifferentialSeq. {(set partials) coefficient}))
  ([coefficient partials & cps]
     (if (odd? (count cps))
       (throw (Exception. "differential-seq requires an even number of terms."))
       (DifferentialSeq.
        (reduce
         #(assoc %1 (set (second %2)) (first %2))
         {(set partials) coefficient}
         (partition 2 cps))))))
  


(defn big-part
  "Returns the part of the differential sequence that is finite,
  i.e. not infinitely small, as a differential sequence. If given a
  non-differential object, returns that object."
  [dseq]
  (if (not= (type dseq) DifferentialSeq) dseq
      (let [m (.terms dseq)
            keys (sort-by count (keys m))
            smallest-var (last (last keys))]
        (DifferentialSeq.
         (reduce
          #(assoc %1 (first %2) (second %2))
          {}
          (remove #((first %) smallest-var) m))))))


(defn small-part
  "Returns the part of the differential sequence that infinitely
  small, as a differential sequence. If given a non-differential
  object, returns that object."
  [dseq]
  (if (not= (type dseq) DifferentialSeq) 0
      (let [m (.terms dseq)
            keys (sort-by count (keys m))
            smallest-var (last (last keys))]
        (DifferentialSeq.
         (reduce
          #(assoc %1 (first %2) (second %2)) {}
          (filter #((first %) smallest-var) m))))))



(defn cartesian-product [set1 set2]
  (reduce concat
          (for [x set1]
            (for [y set2]
              [x y]))))

(defn differential-multiply
  "Multiply two differential sequences. The square of any differential
  variable is zero since differential variables are infinitesimally
  small."
  [dseq1 dseq2]
  (DifferentialSeq.
   (reduce
    (fn [m [[vars1 coeff1] [vars2 coeff2]]]
      (if (not (empty? (clojure.set/intersection vars1 vars2)))
        m
        (assoc m (clojure.set/union vars1 vars2) (* coeff1 coeff2))))
    {}
    (cartesian-product (.terms dseq1) (.terms dseq2)))))



(defmethod * [DifferentialSeq DifferentialSeq]
  [dseq1 dseq2]
   (differential-multiply dseq1 dseq2))

(defmethod + [DifferentialSeq DifferentialSeq]
  [dseq1 dseq2]
  (DifferentialSeq.
   (merge-with + (.terms dseq1) (.terms dseq2))))

(defmethod * [java.lang.Number DifferentialSeq]
  [x dseq]
  (fmap (partial * x) dseq))

(defmethod * [DifferentialSeq java.lang.Number]
  [dseq x]
  (fmap (partial * x) dseq))
;; DERIVATIVES



(defn linear-approximator
  "Returns an operator that linearly approximates the given function."
  ([f df|dx]
     (fn [x]
       (let [big-part (big-part x)
             small-part (small-part x)]
         (+ (fmap f big-part)
            (* (fmap df|dx big-part) small-part)))))
  ([f df|dx df|dy]
     (fn [x y]
       (let [X (big-part x)
             Y (big-part y)
             DX (small-part x)
             DY (small-part y)]
         (+ (f X Y)
            (+ (* DX (df|dx X Y))
               (* DY (df|dy X Y)))))))
  )





(defn D[f]
  (fn[x]
    (f (differential-seq x [] 1 [0]))))

(defn d[f partials]
  (fn [x]
    (get 
     (.terms ((D f)x))
     (set partials)
     0
    )))

(defmethod exp DifferentialSeq [x]
           ((linear-approximator exp exp) x))

(defmethod sin DifferentialSeq
  [x]
  ((linear-approximator sin cos) x))

(defmethod cos DifferentialSeq
  [x]
  ((linear-approximator cos #(- (sin %))) x))
     


Date: 2011-08-11 04:10:36 EDT

Author: Dylan Holmes

Org version 7.6 with Emacs version 23

Validate XHTML 1.0