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))