rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: Building a Classical Mechanics Library in Clojure rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: rlm@2:
rlm@2: rlm@2: rlm@2: rlm@2:
rlm@2:
rlm@2: rlm@2:
rlm@2: rlm@2:

aurellem ☉

rlm@2: rlm@2:
rlm@2: rlm@2:

Building a Classical Mechanics Library in Clojure

rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: rlm@2:
rlm@2:

Table of Contents

rlm@2:
rlm@2: rlm@2:
rlm@2:
rlm@2: rlm@2:
rlm@2:

1 Useful data types

rlm@2:
rlm@2: rlm@2: rlm@2: rlm@2:
rlm@2: rlm@2:
rlm@2:

1.1 Complex numbers

rlm@2:
rlm@2: rlm@2: rlm@2:
rlm@2: rlm@2:
rlm@2: rlm@2:
rlm@2:

1.2 Power series

rlm@2:
rlm@2: rlm@2: rlm@2:
rlm@2: rlm@2:
rlm@2: rlm@2:
rlm@2:

1.3 Tuples and tensors

rlm@2:
rlm@2: rlm@2: rlm@2: rlm@2:
rlm@2: rlm@2:
rlm@2:

1.3.1 Tuples are “sequences with spin”

rlm@2:
rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: 
(in-ns 'sicm.utils)
rlm@2: 
rlm@2: ;; Let some objects have spin
rlm@2: 
rlm@2: (defprotocol Spinning
rlm@2:   (up? [this])
rlm@2:   (down? [this]))
rlm@2: 
rlm@2: (defn spin
rlm@2:   "Returns the spin of the Spinning s, either :up or :down"
rlm@2:   [#^Spinning s]
rlm@2:   (cond (up? s) :up (down? s) :down))
rlm@2: 
rlm@2: 
rlm@2: ;; DEFINITION: A tuple is a sequence with spin
rlm@2: 
rlm@2: (deftype Tuple
rlm@2:   [spin coll]
rlm@2: 
rlm@2:   clojure.lang.Seqable
rlm@2:   (seq [this] (seq (.coll this)))
rlm@2: 
rlm@2:   clojure.lang.Counted
rlm@2:   (count [this] (count (.coll this)))
rlm@2: 
rlm@2:   Spinning
rlm@2:   (up? [this] (= ::up (.spin this)))
rlm@2:   (down? [this] (= ::down (.spin this))))
rlm@2: 
rlm@2: (defmethod print-method Tuple
rlm@2:   [o w]
rlm@2:   (print-simple
rlm@2:    (if (up? o)
rlm@2:      (str "u" (.coll o))
rlm@2:      (str "d" (vec(.coll o))))
rlm@2:    w))
rlm@2: 
rlm@2: 
rlm@2: ;; CONSTRUCTORS
rlm@2: 
rlm@2: (defn up
rlm@2:   "Create a new up-tuple containing the contents of coll."
rlm@2:   [coll]
rlm@2:   (Tuple. ::up coll))       
rlm@2: 
rlm@2: (defn down
rlm@2:   "Create a new down-tuple containing the contents of coll."
rlm@2:   [coll]
rlm@2:   (Tuple. ::down coll))
rlm@2:
rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: rlm@2:
rlm@2: rlm@2:
rlm@2: rlm@2:
rlm@2:

1.3.2 Matrices

rlm@2:
rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: 
(in-ns 'sicm.utils)
rlm@2: (require 'incanter.core) ;; use incanter's fast matrices
rlm@2: 
rlm@2: 
rlm@2: (defn all-equal? [coll]
rlm@2:   (if (empty? (rest coll)) true
rlm@2:       (and (= (first coll) (second coll))
rlm@2:            (recur (rest coll)))))
rlm@2: 
rlm@2: 
rlm@2: (defprotocol Matrix
rlm@2:   (rows [matrix])
rlm@2:   (cols [matrix])
rlm@2:   (diagonal [matrix])
rlm@2:   (trace [matrix])
rlm@2:   (determinant [matrix])
rlm@2:   (transpose [matrix])
rlm@2:   (conjugate [matrix])
rlm@2: )
rlm@2: 
rlm@2: (extend-protocol Matrix incanter.Matrix
rlm@2:   (rows [rs] (map down (apply map vector (apply map vector rs))))
rlm@2:   (cols [rs] (map up (apply map vector rs)))
rlm@2:   (diagonal [matrix] (incanter.core/diag matrix) )
rlm@2:   (determinant [matrix] (incanter.core/det matrix))
rlm@2:   (trace [matrix] (incanter.core/trace matrix))
rlm@2:   (transpose [matrix] (incanter.core/trans matrix)))
rlm@2: 
rlm@2: (defn count-rows [matrix]
rlm@2:   ((comp count rows) matrix))
rlm@2: 
rlm@2: (defn count-cols [matrix]
rlm@2:   ((comp count cols) matrix))
rlm@2: 
rlm@2: (defn square? [matrix]
rlm@2:   (= (count-rows matrix) (count-cols matrix)))
rlm@2: 
rlm@2: (defn identity-matrix
rlm@2:   "Define a square matrix of size n-by-n with 1s along the diagonal and
rlm@2:   0s everywhere else."
rlm@2:   [n]
rlm@2:   (incanter.core/identity-matrix n))
rlm@2: 
rlm@2: 
rlm@2: (defn matrix-by-rows
rlm@2:   "Define a matrix by giving its rows."
rlm@2:   [& rows]
rlm@2:   (if
rlm@2:    (not (all-equal? (map count rows)))
rlm@2:    (throw (Exception. "All rows in a matrix must have the same number of elements."))
rlm@2:    (incanter.core/matrix (vec rows))))
rlm@2: 
rlm@2: (defn matrix-by-cols
rlm@2:   "Define a matrix by giving its columns"
rlm@2:   [& cols]
rlm@2:   (if (not (all-equal? (map count cols)))
rlm@2:    (throw (Exception. "All columns in a matrix must have the same number of elements."))
rlm@2:    (incanter.core/matrix (vec (apply map vector cols)))))
rlm@2: 
rlm@2: (defn identity-matrix
rlm@2:   "Define a square matrix of size n-by-n with 1s along the diagonal and
rlm@2:   0s everywhere else."
rlm@2:   [n]
rlm@2:   (incanter.core/identity-matrix n))
rlm@2: 
rlm@2:
rlm@2: rlm@2: rlm@2: rlm@2: rlm@2:
rlm@2:
rlm@2: rlm@2:
rlm@2: rlm@2:
rlm@2:

1.4 Generic Operations

rlm@2:
rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: 
(in-ns 'sicm.utils)
rlm@2: (use 'clojure.contrib.generic.arithmetic
rlm@2:      'clojure.contrib.generic.collection
rlm@2:      'clojure.contrib.generic.functor
rlm@2:      'clojure.contrib.generic.math-functions)
rlm@2: 
rlm@2: (defn numbers?
rlm@2:   "Returns true if all arguments are numbers, else false."
rlm@2:   [& xs]
rlm@2:   (every? number? xs))
rlm@2: 
rlm@2: (defn tuple-surgery
rlm@2:   "Applies the function f to the items of tuple and the additional
rlm@2:   arguments, if any. Returns a Tuple of the same type as tuple."
rlm@2:   [tuple f & xs]
rlm@2:   ((if (up? tuple) up down)
rlm@2:    (apply f (seq tuple) xs)))
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: ;;; CONTRACTION collapses two compatible tuples into a number.
rlm@2: 
rlm@2: (defn contractible?
rlm@2:   "Returns true if the tuples a and b are compatible for contraction,
rlm@2:   else false. Tuples are compatible if they have the same number of
rlm@2:   components, they have opposite spins, and their elements are
rlm@2:   pairwise-compatible."
rlm@2:   [a b]
rlm@2:   (and
rlm@2:    (isa? (type a) Tuple)
rlm@2:    (isa? (type b) Tuple)
rlm@2:    (= (count a) (count b))
rlm@2:    (not= (spin a) (spin b))
rlm@2:    
rlm@2:    (not-any? false?
rlm@2:              (map #(or
rlm@2:                     (numbers? %1 %2)
rlm@2:                     (contractible? %1 %2))
rlm@2:                   a b))))
rlm@2: 
rlm@2: (defn contract
rlm@2:   "Contracts two tuples, returning the sum of the
rlm@2:   products of the corresponding items. Contraction is recursive on
rlm@2:   nested tuples."
rlm@2:   [a b]
rlm@2:   (if (not (contractible? a b))
rlm@2:     (throw
rlm@2:      (Exception. "Not compatible for contraction."))
rlm@2:     (reduce +
rlm@2:             (map
rlm@2:              (fn [x y]
rlm@2:                (if (numbers? x y)
rlm@2:                  (* x y)
rlm@2:                  (contract x y)))
rlm@2:              a b))))
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: (defmethod conj Tuple
rlm@2:   [tuple & xs]
rlm@2:   (tuple-surgery tuple #(apply conj % xs)))
rlm@2: 
rlm@2: (defmethod fmap Tuple
rlm@2:   [f tuple]
rlm@2:   (tuple-surgery tuple (partial map f)))
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: ;; TODO: define Scalar, and add it to the hierarchy above Number and Complex
rlm@2: 
rlm@2:                                         
rlm@2: (defmethod * [Tuple Tuple]             ; tuple*tuple
rlm@2:   [a b]
rlm@2:   (if (contractible? a b)
rlm@2:     (contract a b)
rlm@2:     (map (partial * a) b)))
rlm@2: 
rlm@2: 
rlm@2: (defmethod * [java.lang.Number Tuple]  ;; scalar *  tuple
rlm@2:   [a x] (fmap (partial * a) x))
rlm@2: 
rlm@2: (defmethod * [Tuple java.lang.Number]
rlm@2:   [x a] (* a x))
rlm@2: 
rlm@2: (defmethod * [java.lang.Number incanter.Matrix] ;; scalar *  matrix
rlm@2:   [x M] (incanter.core/mult x M))
rlm@2: 
rlm@2: (defmethod * [incanter.Matrix java.lang.Number]
rlm@2:   [M x] (* x M))
rlm@2: 
rlm@2: (defmethod * [incanter.Matrix incanter.Matrix] ;; matrix * matrix
rlm@2:   [M1 M2]
rlm@2:   (incanter.core/mmult M1 M2))
rlm@2: 
rlm@2: (defmethod * [incanter.Matrix Tuple] ;; matrix * tuple
rlm@2:   [M v]
rlm@2:   (if (and (apply numbers? v) (up? v)) 
rlm@2:     (* M (matrix-by-cols v))
rlm@2:     (throw (Exception. "Currently, you can only multiply a matrix by a tuple of *numbers*"))
rlm@2:     ))
rlm@2: 
rlm@2: (defmethod * [Tuple incanter.Matrix] ;; tuple * Matrix
rlm@2:   [v M]
rlm@2:   (if (and (apply numbers? v) (down? v))
rlm@2:     (* (matrix-by-rows v) M)
rlm@2:     (throw (Exception. "Currently, you can only multiply a matrix by a tuple of *numbers*"))
rlm@2:     ))
rlm@2: 
rlm@2: 
rlm@2: (defmethod exp incanter.Matrix
rlm@2:   [M]
rlm@2:   (incanter.core/exp M))
rlm@2: 
rlm@2:
rlm@2: rlm@2: rlm@2: rlm@2: rlm@2:
rlm@2:
rlm@2: rlm@2:
rlm@2: rlm@2:
rlm@2:

2 Operators and Differentiation

rlm@2:
rlm@2: rlm@2: rlm@2: rlm@2: rlm@2:
rlm@2: rlm@2:
rlm@2:

2.1 Operators

rlm@2:
rlm@2: rlm@2: rlm@2:
rlm@2: rlm@2:
rlm@2: rlm@2:
rlm@2:

2.2 Differential Terms and Sequences

rlm@2:
rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: 
(in-ns 'sicm.utils)
rlm@2: (use 'clojure.contrib.seq
rlm@2:      'clojure.contrib.generic.arithmetic
rlm@2:      'clojure.contrib.generic.collection
rlm@2:      'clojure.contrib.generic.math-functions)
rlm@2: 
rlm@2: ;;∂
rlm@2: 
rlm@2: ;; DEFINITION : Differential Term
rlm@2: 
rlm@2: ;; A quantity with infinitesimal components, e.g. x, dxdy, 4dydz. The
rlm@2: ;; coefficient of the quantity is returned by the 'coefficient' method,
rlm@2: ;; while the sequence of differential parameters is returned by the
rlm@2: ;; method 'partials'.
rlm@2: 
rlm@2: ;; Instead of using (potentially ambiguous) letters to denote
rlm@2: ;; differential parameters (dx,dy,dz), we use integers. So, dxdz becomes [0 2].
rlm@2: 
rlm@2: ;; The coefficient can be any arithmetic object; the
rlm@2: ;; partials must be a nonrepeating sorted sequence of nonnegative
rlm@2: ;; integers.
rlm@2: 
rlm@2: (deftype DifferentialTerm [coefficient partials])
rlm@2: 
rlm@2: (defn differential-term
rlm@2:   "Make a differential term given pairs of coefficients and partials."
rlm@2:   [coefficient partials]
rlm@2:   (if (and (coll? partials) (every? #(and (integer? %) (not(neg? %))) partials)) 
rlm@2:     (DifferentialTerm. coefficient (set partials))
rlm@2:     (throw (java.lang.IllegalArgumentException. "Partials must be a collection of integers."))))
rlm@2: 
rlm@2: 
rlm@2: ;; DEFINITION : Differential Sequence
rlm@2: ;; A differential sequence is a sequence of differential terms, all with different partials.
rlm@2: ;; Internally, it is a map from the partials of each term to their coefficients.
rlm@2: 
rlm@2: (deftype DifferentialSeq
rlm@2:   [terms]
rlm@2:   ;;clojure.lang.IPersistentMap
rlm@2:   clojure.lang.Associative
rlm@2:   (assoc [this key val]
rlm@2:     (DifferentialSeq.
rlm@2:      (cons (differential-term val key) terms)))
rlm@2:   (cons [this x]
rlm@2:         (DifferentialSeq. (cons x terms)))
rlm@2:   (containsKey [this key]
rlm@2:                (not(nil? (find-first #(= (.partials %) key) terms))))
rlm@2:   (count [this] (count (.terms this)))
rlm@2:   (empty [this] (DifferentialSeq. []))
rlm@2:   (entryAt [this key]
rlm@2:            ((juxt #(.partials %) #(.coefficient %))
rlm@2:             (find-first #(= (.partials %) key) terms)))
rlm@2:   (seq [this] (seq (.terms this))))
rlm@2: 
rlm@2: 
rlm@2: (defmethod fmap DifferentialSeq
rlm@2:   [f dseq]
rlm@2:   (DifferentialSeq. (fmap f (.terms dseq))))
rlm@2: 
rlm@2: 
rlm@2: ;; BUILDING DIFFERENTIAL OBJECTS
rlm@2: 
rlm@2: (defn differential-seq
rlm@2:     "Define a differential sequence by specifying coefficient/partials
rlm@2: pairs, which are used to create differential terms to populate the
rlm@2: sequence."
rlm@2:   ([coefficient partials]
rlm@2:      (DifferentialSeq. {(set partials) coefficient}))
rlm@2:   ([coefficient partials & cps]
rlm@2:      (if (odd? (count cps))
rlm@2:        (throw (Exception. "differential-seq requires an even number of terms."))
rlm@2:        (DifferentialSeq.
rlm@2:         (reduce
rlm@2:          #(assoc %1 (set (second %2)) (first %2))
rlm@2:          {(set partials) coefficient}
rlm@2:          (partition 2 cps))))))
rlm@2:   
rlm@2: 
rlm@2: 
rlm@2: (defn big-part
rlm@2:   "Returns the part of the differential sequence that is finite,
rlm@2:   i.e. not infinitely small, as a differential sequence. If given a
rlm@2:   non-differential object, returns that object."
rlm@2:   [dseq]
rlm@2:   (if (not= (type dseq) DifferentialSeq) dseq
rlm@2:       (let [m (.terms dseq)
rlm@2:             keys (sort-by count (keys m))
rlm@2:             smallest-var (last (last keys))]
rlm@2:         (DifferentialSeq.
rlm@2:          (reduce
rlm@2:           #(assoc %1 (first %2) (second %2))
rlm@2:           {}
rlm@2:           (remove #((first %) smallest-var) m))))))
rlm@2: 
rlm@2: 
rlm@2: (defn small-part
rlm@2:   "Returns the part of the differential sequence that infinitely
rlm@2:   small, as a differential sequence. If given a non-differential
rlm@2:   object, returns that object."
rlm@2:   [dseq]
rlm@2:   (if (not= (type dseq) DifferentialSeq) 0
rlm@2:       (let [m (.terms dseq)
rlm@2:             keys (sort-by count (keys m))
rlm@2:             smallest-var (last (last keys))]
rlm@2:         (DifferentialSeq.
rlm@2:          (reduce
rlm@2:           #(assoc %1 (first %2) (second %2)) {}
rlm@2:           (filter #((first %) smallest-var) m))))))
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: (defn cartesian-product [set1 set2]
rlm@2:   (reduce concat
rlm@2:           (for [x set1]
rlm@2:             (for [y set2]
rlm@2:               [x y]))))
rlm@2: 
rlm@2: (defn differential-multiply
rlm@2:   "Multiply two differential sequences. The square of any differential
rlm@2:   variable is zero since differential variables are infinitesimally
rlm@2:   small."
rlm@2:   [dseq1 dseq2]
rlm@2:   (DifferentialSeq.
rlm@2:    (reduce
rlm@2:     (fn [m [[vars1 coeff1] [vars2 coeff2]]]
rlm@2:       (if (not (empty? (clojure.set/intersection vars1 vars2)))
rlm@2:         m
rlm@2:         (assoc m (clojure.set/union vars1 vars2) (* coeff1 coeff2))))
rlm@2:     {}
rlm@2:     (cartesian-product (.terms dseq1) (.terms dseq2)))))
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: (defmethod * [DifferentialSeq DifferentialSeq]
rlm@2:   [dseq1 dseq2]
rlm@2:    (differential-multiply dseq1 dseq2))
rlm@2: 
rlm@2: (defmethod + [DifferentialSeq DifferentialSeq]
rlm@2:   [dseq1 dseq2]
rlm@2:   (DifferentialSeq.
rlm@2:    (merge-with + (.terms dseq1) (.terms dseq2))))
rlm@2: 
rlm@2: (defmethod * [java.lang.Number DifferentialSeq]
rlm@2:   [x dseq]
rlm@2:   (fmap (partial * x) dseq))
rlm@2: 
rlm@2: (defmethod * [DifferentialSeq java.lang.Number]
rlm@2:   [dseq x]
rlm@2:   (fmap (partial * x) dseq))
rlm@2: ;; DERIVATIVES
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: (defn linear-approximator
rlm@2:   "Returns an operator that linearly approximates the given function."
rlm@2:   ([f df|dx]
rlm@2:      (fn [x]
rlm@2:        (let [big-part (big-part x)
rlm@2:              small-part (small-part x)]
rlm@2:          (+ (fmap f big-part)
rlm@2:             (* (fmap df|dx big-part) small-part)))))
rlm@2:   ([f df|dx df|dy]
rlm@2:      (fn [x y]
rlm@2:        (let [X (big-part x)
rlm@2:              Y (big-part y)
rlm@2:              DX (small-part x)
rlm@2:              DY (small-part y)]
rlm@2:          (+ (f X Y)
rlm@2:             (+ (* DX (df|dx X Y))
rlm@2:                (* DY (df|dy X Y)))))))
rlm@2:   )
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: (defn D[f]
rlm@2:   (fn[x]
rlm@2:     (f (differential-seq x [] 1 [0]))))
rlm@2: 
rlm@2: (defn d[f partials]
rlm@2:   (fn [x]
rlm@2:     (get 
rlm@2:      (.terms ((D f)x))
rlm@2:      (set partials)
rlm@2:      0
rlm@2:     )))
rlm@2: 
rlm@2: (defmethod exp DifferentialSeq [x]
rlm@2:            ((linear-approximator exp exp) x))
rlm@2: 
rlm@2: (defmethod sin DifferentialSeq
rlm@2:   [x]
rlm@2:   ((linear-approximator sin cos) x))
rlm@2: 
rlm@2: (defmethod cos DifferentialSeq
rlm@2:   [x]
rlm@2:   ((linear-approximator cos #(- (sin %))) x))
rlm@2:      
rlm@2: 
rlm@2: 
rlm@2:
rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: rlm@2: rlm@2:
rlm@2:
rlm@2:
rlm@2:
rlm@2:

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

rlm@2:

Author: Dylan Holmes

rlm@2:

Org version 7.6 with Emacs version 23

rlm@2: Validate XHTML 1.0 rlm@2:
rlm@2:
rlm@2: rlm@2: