Mercurial > dylan

annotate sicm/utils.org @ 5:3ff40c869d1a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
ReTouchups
author Dylan Holmes <ocsenave@gmail.com>
date Thu, 27 Oct 2011 22:38:23 -0500
parents b4de894a1e2e
children
rev   line source
rlm@2 1 #+TITLE: Building a Classical Mechanics Library in Clojure
rlm@2 2 #+author: Dylan Holmes
rlm@2 3 #+EMAIL: rlm@mit.edu
rlm@2 4 #+MATHJAX: align:"left" mathml:t path:"../MathJax/MathJax.js"
rlm@2 5 #+STYLE: <link rel="stylesheet" type="text/css" href="../css/aurellem.css" />
rlm@2 6 #+OPTIONS: H:3 num:t toc:t \n:nil @:t ::t |:t ^:t -:t f:t *:t <:t
rlm@2 7 #+SETUPFILE: ../templates/level-0.org
rlm@2 8 #+INCLUDE: ../templates/level-0.org
rlm@2 9 #+BABEL: :noweb yes :results silent
rlm@2 10
rlm@2 11 * Useful data types
rlm@2 12
rlm@2 13 ** Complex numbers
rlm@2 14
rlm@2 15 ** Power series
rlm@2 16
rlm@2 17 ** Tuples and tensors
rlm@2 18
rlm@2 19 *** Tuples are \ldquo{}sequences with spin\rdquo{}
rlm@2 20
rlm@2 21 #+srcname: tuples
rlm@2 22 #+begin_src clojure
rlm@2 23 (in-ns 'sicm.utils)
rlm@2 24
rlm@2 25 ;; Let some objects have spin
rlm@2 26
rlm@2 27 (defprotocol Spinning
rlm@2 28 (up? [this])
rlm@2 29 (down? [this]))
rlm@2 30
rlm@2 31 (defn spin
rlm@2 32 "Returns the spin of the Spinning s, either :up or :down"
rlm@2 33 [#^Spinning s]
rlm@2 34 (cond (up? s) :up (down? s) :down))
rlm@2 35
rlm@2 36
rlm@2 37 ;; DEFINITION: A tuple is a sequence with spin
rlm@2 38
rlm@2 39 (deftype Tuple
rlm@2 40 [spin coll]
rlm@2 41
rlm@2 42 clojure.lang.Seqable
rlm@2 43 (seq [this] (seq (.coll this)))
rlm@2 44
rlm@2 45 clojure.lang.Counted
rlm@2 46 (count [this] (count (.coll this)))
rlm@2 47
rlm@2 48 Spinning
rlm@2 49 (up? [this] (= ::up (.spin this)))
rlm@2 50 (down? [this] (= ::down (.spin this))))
rlm@2 51
rlm@2 52 (defmethod print-method Tuple
rlm@2 53 [o w]
rlm@2 54 (print-simple
rlm@2 55 (if (up? o)
rlm@2 56 (str "u" (.coll o))
rlm@2 57 (str "d" (vec(.coll o))))
rlm@2 58 w))
rlm@2 59
rlm@2 60 (def tuple? #(= (type %) Tuple))
rlm@2 61
rlm@2 62 ;; CONSTRUCTORS
rlm@2 63
rlm@2 64 (defn up
rlm@2 65 "Create a new up-tuple containing the contents of coll."
rlm@2 66 [coll]
rlm@2 67 (Tuple. ::up coll))
rlm@2 68
rlm@2 69 (defn down
rlm@2 70 "Create a new down-tuple containing the contents of coll."
rlm@2 71 [coll]
rlm@2 72 (Tuple. ::down coll))
rlm@2 73
rlm@2 74 (defn same-spin
rlm@2 75 "Creates a tuple which has the same spin as tuple and which contains
rlm@2 76 the contents of coll."
rlm@2 77 [tuple coll]
rlm@2 78 (if (up? tuple)
rlm@2 79 (up coll)
rlm@2 80 (down coll)))
rlm@2 81
rlm@2 82 (defn opposite-spin
rlm@2 83 "Create a tuple which has opposite spin to tuple and which contains
rlm@2 84 the contents of coll."
rlm@2 85 [tuple coll]
rlm@2 86 (if (up? tuple)
rlm@2 87 (down coll)
rlm@2 88 (up coll)))
rlm@2 89 #+end_src
rlm@2 90
rlm@2 91
rlm@2 92
rlm@2 93 *** Matrices
rlm@2 94 #+srcname:matrices
rlm@2 95 #+begin_src clojure
rlm@2 96 (in-ns 'sicm.utils)
rlm@2 97 (require 'incanter.core) ;; use incanter's fast matrices
rlm@2 98
rlm@2 99
rlm@2 100 (defn all-equal? [coll]
rlm@2 101 (if (empty? (rest coll)) true
rlm@2 102 (and (= (first coll) (second coll))
rlm@2 103 (recur (rest coll)))))
rlm@2 104
rlm@2 105
rlm@2 106 (defprotocol Matrix
rlm@2 107 (rows [matrix])
rlm@2 108 (cols [matrix])
rlm@2 109 (diagonal [matrix])
rlm@2 110 (trace [matrix])
rlm@2 111 (determinant [matrix])
rlm@2 112 (transpose [matrix])
rlm@2 113 (conjugate [matrix])
rlm@2 114 )
rlm@2 115
rlm@2 116 (extend-protocol Matrix incanter.Matrix
rlm@2 117 (rows [rs] (map down (apply map vector (apply map vector rs))))
rlm@2 118 (cols [rs] (map up (apply map vector rs)))
rlm@2 119 (diagonal [matrix] (incanter.core/diag matrix) )
rlm@2 120 (determinant [matrix] (incanter.core/det matrix))
rlm@2 121 (trace [matrix] (incanter.core/trace matrix))
rlm@2 122 (transpose [matrix] (incanter.core/trans matrix)))
rlm@2 123
rlm@2 124 (defn count-rows [matrix]
rlm@2 125 ((comp count rows) matrix))
rlm@2 126
rlm@2 127 (defn count-cols [matrix]
rlm@2 128 ((comp count cols) matrix))
rlm@2 129
rlm@2 130 (defn square? [matrix]
rlm@2 131 (= (count-rows matrix) (count-cols matrix)))
rlm@2 132
rlm@2 133 (defn identity-matrix
rlm@2 134 "Define a square matrix of size n-by-n with 1s along the diagonal and
rlm@2 135 0s everywhere else."
rlm@2 136 [n]
rlm@2 137 (incanter.core/identity-matrix n))
rlm@2 138
rlm@2 139
rlm@2 140 (defn matrix-by-rows
rlm@2 141 "Define a matrix by giving its rows."
rlm@2 142 [& rows]
rlm@2 143 (if
rlm@2 144 (not (all-equal? (map count rows)))
rlm@2 145 (throw (Exception. "All rows in a matrix must have the same number of elements."))
rlm@2 146 (incanter.core/matrix (vec rows))))
rlm@2 147
rlm@2 148 (defn matrix-by-cols
rlm@2 149 "Define a matrix by giving its columns"
rlm@2 150 [& cols]
rlm@2 151 (if (not (all-equal? (map count cols)))
rlm@2 152 (throw (Exception. "All columns in a matrix must have the same number of elements."))
rlm@2 153 (incanter.core/matrix (vec (apply map vector cols)))))
rlm@2 154
rlm@2 155 (defn identity-matrix
rlm@2 156 "Define a square matrix of size n-by-n with 1s along the diagonal and
rlm@2 157 0s everywhere else."
rlm@2 158 [n]
rlm@2 159 (incanter.core/identity-matrix n))
rlm@2 160
rlm@2 161 #+end_src
rlm@2 162
rlm@2 163 ** Generic Operations
rlm@2 164 #+srcname:arith-tuple
rlm@2 165 #+begin_src clojure
rlm@2 166 (in-ns 'sicm.utils)
rlm@2 167 (use 'clojure.contrib.generic.arithmetic
rlm@2 168 'clojure.contrib.generic.collection
rlm@2 169 'clojure.contrib.generic.functor
rlm@2 170 'clojure.contrib.generic.math-functions)
rlm@2 171
rlm@2 172 (defn numbers?
rlm@2 173 "Returns true if all arguments are numbers, else false."
rlm@2 174 [& xs]
rlm@2 175 (every? number? xs))
rlm@2 176
rlm@2 177 (defn tuple-surgery
rlm@2 178 "Applies the function f to the items of tuple and the additional
rlm@2 179 arguments, if any. Returns a Tuple of the same type as tuple."
rlm@2 180 [tuple f & xs]
rlm@2 181 ((if (up? tuple) up down)
rlm@2 182 (apply f (seq tuple) xs)))
rlm@2 183
rlm@2 184
rlm@2 185
rlm@2 186 ;;; CONTRACTION collapses two compatible tuples into a number.
rlm@2 187
rlm@2 188 (defn contractible?
rlm@2 189 "Returns true if the tuples a and b are compatible for contraction,
rlm@2 190 else false. Tuples are compatible if they have the same number of
rlm@2 191 components, they have opposite spins, and their elements are
rlm@2 192 pairwise-compatible."
rlm@2 193 [a b]
rlm@2 194 (and
rlm@2 195 (isa? (type a) Tuple)
rlm@2 196 (isa? (type b) Tuple)
rlm@2 197 (= (count a) (count b))
rlm@2 198 (not= (spin a) (spin b))
rlm@2 199
rlm@2 200 (not-any? false?
rlm@2 201 (map #(or
rlm@2 202 (numbers? %1 %2)
rlm@2 203 (contractible? %1 %2))
rlm@2 204 a b))))
rlm@2 205
rlm@2 206 (defn contract
rlm@2 207 "Contracts two tuples, returning the sum of the
rlm@2 208 products of the corresponding items. Contraction is recursive on
rlm@2 209 nested tuples."
rlm@2 210 [a b]
rlm@2 211 (if (not (contractible? a b))
rlm@2 212 (throw
rlm@2 213 (Exception. "Not compatible for contraction."))
rlm@2 214 (reduce +
rlm@2 215 (map
rlm@2 216 (fn [x y]
rlm@2 217 (if (numbers? x y)
rlm@2 218 (* x y)
rlm@2 219 (contract x y)))
rlm@2 220 a b))))
rlm@2 221
rlm@2 222
rlm@2 223
rlm@2 224
rlm@2 225
rlm@2 226 (defmethod conj Tuple
rlm@2 227 [tuple & xs]
rlm@2 228 (tuple-surgery tuple #(apply conj % xs)))
rlm@2 229
rlm@2 230 (defmethod fmap Tuple
rlm@2 231 [f tuple]
rlm@2 232 (tuple-surgery tuple (partial map f)))
rlm@2 233
rlm@2 234
rlm@2 235
rlm@2 236 ;; TODO: define Scalar, and add it to the hierarchy above Number and Complex
rlm@2 237
rlm@2 238
rlm@2 239 (defmethod * [Tuple Tuple] ; tuple*tuple
rlm@2 240 [a b]
rlm@2 241 (if (contractible? a b)
rlm@2 242 (contract a b)
rlm@2 243 (map (partial * a) b)))
rlm@2 244
rlm@2 245
rlm@2 246 (defmethod * [java.lang.Number Tuple] ;; scalar * tuple
rlm@2 247 [a x] (fmap (partial * a) x))
rlm@2 248
rlm@2 249 (defmethod * [Tuple java.lang.Number]
rlm@2 250 [x a] (* a x))
rlm@2 251
rlm@2 252 (defmethod * [java.lang.Number incanter.Matrix] ;; scalar * matrix
rlm@2 253 [x M] (incanter.core/mult x M))
rlm@2 254
rlm@2 255 (defmethod * [incanter.Matrix java.lang.Number]
rlm@2 256 [M x] (* x M))
rlm@2 257
rlm@2 258 (defmethod * [incanter.Matrix incanter.Matrix] ;; matrix * matrix
rlm@2 259 [M1 M2]
rlm@2 260 (incanter.core/mmult M1 M2))
rlm@2 261
rlm@2 262 (defmethod * [incanter.Matrix Tuple] ;; matrix * tuple
rlm@2 263 [M v]
rlm@2 264 (if (and (apply numbers? v) (up? v))
rlm@2 265 (* M (matrix-by-cols v))
rlm@2 266 (throw (Exception. "Currently, you can only multiply a matrix by a tuple of *numbers*"))
rlm@2 267 ))
rlm@2 268
rlm@2 269 (defmethod * [Tuple incanter.Matrix] ;; tuple * Matrix
rlm@2 270 [v M]
rlm@2 271 (if (and (apply numbers? v) (down? v))
rlm@2 272 (* (matrix-by-rows v) M)
rlm@2 273 (throw (Exception. "Currently, you can only multiply a matrix by a tuple of *numbers*"))
rlm@2 274 ))
rlm@2 275
rlm@2 276
rlm@2 277 (defmethod exp incanter.Matrix
rlm@2 278 [M]
rlm@2 279 (incanter.core/exp M))
rlm@2 280
rlm@2 281 #+end_src
rlm@2 282
rlm@2 283 * Operators and Differentiation
rlm@2 284 ** Operators
rlm@2 285 #+scrname: operators
rlm@2 286 #+begin_src clojure
rlm@2 287 (in-ns 'sicm.utils)
rlm@2 288 (use 'clojure.contrib.seq
rlm@2 289 'clojure.contrib.generic.arithmetic
rlm@2 290 'clojure.contrib.generic.collection
rlm@2 291 'clojure.contrib.generic.math-functions)
rlm@2 292
rlm@2 293 (defmethod + [clojure.lang.IFn clojure.lang.IFn]
rlm@2 294 [f g]
rlm@2 295 (fn [& args]
rlm@2 296 (+ (apply f args) (apply g args))))
rlm@2 297
rlm@2 298 (defmethod * [clojure.lang.IFn clojure.lang.IFn]
rlm@2 299 [f g]
rlm@2 300 (fn [& args]
rlm@2 301 (* (apply f args) (apply g args))))
rlm@2 302
rlm@2 303 (defmethod / [clojure.lang.IFn java.lang.Number]
rlm@2 304 [f x]
rlm@2 305 (fn [& args]
rlm@2 306 (/ (apply f args) x)))
rlm@2 307
rlm@2 308
rlm@2 309 (defmethod - [clojure.lang.IFn]
rlm@2 310 [f]
rlm@2 311 (fn [& args]
rlm@2 312 (- (apply f args))))
rlm@2 313
rlm@2 314 (defmethod - [clojure.lang.IFn clojure.lang.IFn]
rlm@2 315 [f g]
rlm@2 316 (fn [& args]
rlm@2 317 (- (apply f args) (apply g args))))
rlm@2 318
rlm@2 319 (defmethod pow [clojure.lang.IFn java.lang.Number]
rlm@2 320 [f x]
rlm@2 321 (fn [& args]
rlm@2 322 (pow (apply f args) x)))
rlm@2 323
rlm@2 324
rlm@2 325 (defmethod + [java.lang.Number clojure.lang.IFn]
rlm@2 326 [x f]
rlm@2 327 (fn [& args]
rlm@2 328 (+ x (apply f args))))
rlm@2 329
rlm@2 330 (defmethod * [java.lang.Number clojure.lang.IFn]
rlm@2 331 [x f]
rlm@2 332 (fn [& args]
rlm@2 333 (* x (apply f args))))
rlm@2 334
rlm@2 335 (defmethod * [clojure.lang.IFn java.lang.Number]
rlm@2 336 [f x]
rlm@2 337 (* x f))
rlm@2 338 (defmethod + [clojure.lang.IFn java.lang.Number]
rlm@2 339 [f x]
rlm@2 340 (+ x f))
rlm@2 341
rlm@2 342 #+end_src
rlm@2 343
rlm@2 344 ** Differential Terms and Sequences
rlm@2 345 #+srcname: differential
rlm@2 346 #+begin_src clojure
rlm@2 347 (in-ns 'sicm.utils)
rlm@2 348 (use 'clojure.contrib.seq
rlm@2 349 'clojure.contrib.generic.arithmetic
rlm@2 350 'clojure.contrib.generic.collection
rlm@2 351 'clojure.contrib.generic.math-functions)
rlm@2 352
rlm@2 353 ;;∂
rlm@2 354
rlm@2 355 ;; DEFINITION : Differential Term
rlm@2 356
rlm@2 357 ;; A quantity with infinitesimal components, e.g. x, dxdy, 4dydz. The
rlm@2 358 ;; coefficient of the quantity is returned by the 'coefficient' method,
rlm@2 359 ;; while the sequence of differential parameters is returned by the
rlm@2 360 ;; method 'partials'.
rlm@2 361
rlm@2 362 ;; Instead of using (potentially ambiguous) letters to denote
rlm@2 363 ;; differential parameters (dx,dy,dz), we use integers. So, dxdz becomes [0 2].
rlm@2 364
rlm@2 365 ;; The coefficient can be any arithmetic object; the
rlm@2 366 ;; partials must be a nonrepeating sorted sequence of nonnegative
rlm@2 367 ;; integers.
rlm@2 368
rlm@2 369 ;; (deftype DifferentialTerm [coefficient partials])
rlm@2 370
rlm@2 371 ;; (defn differential-term
rlm@2 372 ;; "Make a differential term from a coefficient and list of partials."
rlm@2 373 ;; [coefficient partials]
rlm@2 374 ;; (if (and (coll? partials) (every? #(and (integer? %) (not(neg? %))) partials))
rlm@2 375 ;; (DifferentialTerm. coefficient (set partials))
rlm@2 376 ;; (throw (java.lang.IllegalArgumentException. "Partials must be a collection of integers."))))
rlm@2 377
rlm@2 378
rlm@2 379 ;; DEFINITION : Differential Sequence
rlm@2 380 ;; A differential sequence is a sequence of differential terms, all with different partials.
rlm@2 381 ;; Internally, it is a map from the partials of each term to their coefficients.
rlm@2 382
rlm@2 383 (deftype DifferentialSeq
rlm@2 384 [terms]
rlm@2 385 ;;clojure.lang.IPersistentMap
rlm@2 386 clojure.lang.Associative
rlm@2 387 (assoc [this key val]
rlm@2 388 (DifferentialSeq.
rlm@2 389 (cons (differential-term val key) terms)))
rlm@2 390 (cons [this x]
rlm@2 391 (DifferentialSeq. (cons x terms)))
rlm@2 392 (containsKey [this key]
rlm@2 393 (not(nil? (find-first #(= (.partials %) key) terms))))
rlm@2 394 (count [this] (count (.terms this)))
rlm@2 395 (empty [this] (DifferentialSeq. []))
rlm@2 396 (entryAt [this key]
rlm@2 397 ((juxt #(.partials %) #(.coefficient %))
rlm@2 398 (find-first #(= (.partials %) key) terms)))
rlm@2 399 (seq [this] (seq (.terms this))))
rlm@2 400
rlm@2 401 (def differential? #(= (type %) DifferentialSeq))
rlm@2 402
rlm@2 403 (defn zeroth-order?
rlm@2 404 "Returns true if the differential sequence has at most a constant term."
rlm@2 405 [dseq]
rlm@2 406 (and
rlm@2 407 (differential? dseq)
rlm@2 408 (every?
rlm@2 409 #(= #{} %)
rlm@2 410 (keys (.terms dseq)))))
rlm@2 411
rlm@2 412 (defmethod fmap DifferentialSeq
rlm@2 413 [f dseq]
rlm@2 414 (DifferentialSeq.
rlm@2 415 (fmap f (.terms dseq))))
rlm@2 416
rlm@2 417
rlm@2 418
rlm@2 419
rlm@2 420 ;; BUILDING DIFFERENTIAL OBJECTS
rlm@2 421
rlm@2 422 (defn differential-seq
rlm@2 423 "Define a differential sequence by specifying an alternating
rlm@2 424 sequence of coefficients and lists of partials."
rlm@2 425 ([coefficient partials]
rlm@2 426 (DifferentialSeq. {(set partials) coefficient}))
rlm@2 427 ([coefficient partials & cps]
rlm@2 428 (if (odd? (count cps))
rlm@2 429 (throw (Exception. "differential-seq requires an even number of terms."))
rlm@2 430 (DifferentialSeq.
rlm@2 431 (reduce
rlm@2 432 #(assoc %1 (set (second %2)) (first %2))
rlm@2 433 {(set partials) coefficient}
rlm@2 434 (partition 2 cps))))))
rlm@2 435
rlm@2 436
rlm@2 437
rlm@2 438 (defn big-part
rlm@2 439 "Returns the part of the differential sequence that is finite,
rlm@2 440 i.e. not infinitely small. If the sequence is zeroth-order, returns
rlm@2 441 the coefficient of the zeroth-order term instead. "
rlm@2 442 [dseq]
rlm@2 443 (if (zeroth-order? dseq) (get (.terms dseq) #{})
rlm@2 444 (let [m (.terms dseq)
rlm@2 445 keys (sort-by count (keys m))
rlm@2 446 smallest-var (last (last keys))]
rlm@2 447 (DifferentialSeq.
rlm@2 448 (reduce
rlm@2 449 #(assoc %1 (first %2) (second %2))
rlm@2 450 {}
rlm@2 451 (remove #((first %) smallest-var) m))))))
rlm@2 452
rlm@2 453
rlm@2 454 (defn small-part
rlm@2 455 "Returns the part of the differential sequence that infinitely
rlm@2 456 small. If the sequence is zeroth-order, returns zero."
rlm@2 457 [dseq]
rlm@2 458 (if (zeroth-order? dseq) 0
rlm@2 459 (let [m (.terms dseq)
rlm@2 460 keys (sort-by count (keys m))
rlm@2 461 smallest-var (last (last keys))]
rlm@2 462 (DifferentialSeq.
rlm@2 463 (reduce
rlm@2 464 #(assoc %1 (first %2) (second %2)) {}
rlm@2 465 (filter #((first %) smallest-var) m))))))
rlm@2 466
rlm@2 467
rlm@2 468
rlm@2 469 (defn cartesian-product [set1 set2]
rlm@2 470 (reduce concat
rlm@2 471 (for [x set1]
rlm@2 472 (for [y set2]
rlm@2 473 [x y]))))
rlm@2 474
rlm@2 475 (defn nth-subset [n]
rlm@2 476 (if (zero? n) []
rlm@2 477 (let [lg2 #(/ (log %) (log 2))
rlm@2 478 k (int(java.lang.Math/floor (lg2 n)))
rlm@2 479 ]
rlm@2 480 (cons k
rlm@2 481 (nth-subset (- n (pow 2 k)))))))
rlm@2 482
rlm@2 483 (def all-partials
rlm@2 484 (lazy-seq (map nth-subset (range))))
rlm@2 485
rlm@2 486
rlm@2 487 (defn differential-multiply
rlm@2 488 "Multiply two differential sequences. The square of any differential
rlm@2 489 variable is zero since differential variables are infinitesimally
rlm@2 490 small."
rlm@2 491 [dseq1 dseq2]
rlm@2 492 (DifferentialSeq.
rlm@2 493 (reduce
rlm@2 494 (fn [m [[vars1 coeff1] [vars2 coeff2]]]
rlm@2 495 (if (not (empty? (clojure.set/intersection vars1 vars2)))
rlm@2 496 m
rlm@2 497 (assoc m (clojure.set/union vars1 vars2) (* coeff1 coeff2))))
rlm@2 498 {}
rlm@2 499 (cartesian-product (.terms dseq1) (.terms dseq2)))))
rlm@2 500
rlm@2 501
rlm@2 502
rlm@2 503 (defmethod * [DifferentialSeq DifferentialSeq]
rlm@2 504 [dseq1 dseq2]
rlm@2 505 (differential-multiply dseq1 dseq2))
rlm@2 506
rlm@2 507 (defmethod + [DifferentialSeq DifferentialSeq]
rlm@2 508 [dseq1 dseq2]
rlm@2 509 (DifferentialSeq.
rlm@2 510 (merge-with + (.terms dseq1) (.terms dseq2))))
rlm@2 511
rlm@2 512 (defmethod * [java.lang.Number DifferentialSeq]
rlm@2 513 [x dseq]
rlm@2 514 (fmap (partial * x) dseq))
rlm@2 515
rlm@2 516 (defmethod * [DifferentialSeq java.lang.Number]
rlm@2 517 [dseq x]
rlm@2 518 (fmap (partial * x) dseq))
rlm@2 519
rlm@2 520 (defmethod + [java.lang.Number DifferentialSeq]
rlm@2 521 [x dseq]
rlm@2 522 (+ (differential-seq x []) dseq))
rlm@2 523 (defmethod + [DifferentialSeq java.lang.Number]
rlm@2 524 [dseq x]
rlm@2 525 (+ dseq (differential-seq x [])))
rlm@2 526
rlm@2 527 (defmethod - DifferentialSeq
rlm@2 528 [x]
rlm@2 529 (fmap - x))
rlm@2 530
rlm@2 531
rlm@2 532 ;; DERIVATIVES
rlm@2 533
rlm@2 534
rlm@2 535
rlm@2 536 (defn linear-approximator
rlm@2 537 "Returns an operator that linearly approximates the given function."
rlm@2 538 ([f df|dx]
rlm@2 539 (fn [x]
rlm@2 540 (let [big-part (big-part x)
rlm@2 541 small-part (small-part x)]
rlm@2 542 ;; f(x+dx) ~= f(x) + f'(x)dx
rlm@2 543 (+ (f big-part)
rlm@2 544 (* (df|dx big-part) small-part)
rlm@2 545 ))))
rlm@2 546
rlm@2 547 ([f df|dx df|dy]
rlm@2 548 (fn [x y]
rlm@2 549 (let [X (big-part x)
rlm@2 550 Y (big-part y)
rlm@2 551 DX (small-part x)
rlm@2 552 DY (small-part y)]
rlm@2 553 (+ (f X Y)
rlm@2 554 (* DX (f df|dx X Y))
rlm@2 555 (* DY (f df|dy X Y)))))))
rlm@2 556
rlm@2 557
rlm@2 558
rlm@2 559
rlm@2 560
rlm@2 561 (defn D[f]
rlm@2 562 (fn[x] (f (+ x (differential-seq 1 [0] 1 [1] 1 [2])))))
rlm@2 563
rlm@2 564 (defn d[partials f]
rlm@2 565 (fn [x]
rlm@2 566 (get
rlm@2 567 (.terms ((D f)x))
rlm@2 568 (set partials)
rlm@2 569 0
rlm@2 570 )))
rlm@2 571
rlm@2 572 (defmethod exp DifferentialSeq [x]
rlm@2 573 ((linear-approximator exp exp) x))
rlm@2 574
rlm@2 575 (defmethod sin DifferentialSeq
rlm@2 576 [x]
rlm@2 577 ((linear-approximator sin cos) x))
rlm@2 578
rlm@2 579 (defmethod cos DifferentialSeq
rlm@2 580 [x]
rlm@2 581 ((linear-approximator cos #(- (sin %))) x))
rlm@2 582
rlm@2 583 (defmethod log DifferentialSeq
rlm@2 584 [x]
rlm@2 585 ((linear-approximator log (fn [x] (/ x)) ) x))
rlm@2 586
rlm@2 587 (defmethod / [DifferentialSeq DifferentialSeq]
rlm@2 588 [x y]
rlm@2 589 ((linear-approximator /
rlm@2 590 (fn [x y] (/ 1 y))
rlm@2 591 (fn [x y] (- (/ x (* y y)))))
rlm@2 592 x y))
rlm@2 593
rlm@2 594 #+end_src
rlm@2 595
rlm@2 596
rlm@2 597
rlm@2 598 * Derivatives Revisited
rlm@2 599 #+begin_src clojure
rlm@2 600 (in-ns 'sicm.utils)
rlm@2 601 (use 'clojure.contrib.seq
rlm@2 602 'clojure.contrib.generic.arithmetic
rlm@2 603 'clojure.contrib.generic.collection
rlm@2 604 'clojure.contrib.generic.math-functions)
rlm@2 605
rlm@2 606 (defn replace-in
rlm@2 607 "Replaces the nth item in coll with the given item. If n is larger
rlm@2 608 than the size of coll, adds n to the end of the collection."
rlm@2 609 [coll n item]
rlm@2 610 (concat
rlm@2 611 (take n coll)
rlm@2 612 [item]
rlm@2 613 (drop (inc n) coll)))
rlm@2 614
rlm@2 615 (defn euclidean-structure [f partials]
rlm@2 616 (fn sd[g v]
rlm@2 617 (cond
rlm@2 618 (tuple? v)
rlm@2 619 (opposite-spin
rlm@2 620 v
rlm@2 621 (map
rlm@2 622 (fn [n]
rlm@2 623 (sd (fn [xn]
rlm@2 624 (g
rlm@2 625 (same-spin v
rlm@2 626 (replace-in v n xn))))
rlm@2 627 (nth v n)))
rlm@2 628 (range (count v)))))))
rlm@2 629
rlm@2 630
rlm@2 631
rlm@2 632
rlm@2 633 #+end_src
rlm@2 634
rlm@2 635
rlm@2 636 * Symbolic Quantities
rlm@2 637
rlm@2 638 #+srcname: symbolic
rlm@2 639 #+begin_src clojure
rlm@2 640 (in-ns 'sicm.utils)
rlm@2 641 (use 'clojure.contrib.generic.arithmetic
rlm@2 642 'clojure.contrib.generic.math-functions)
rlm@2 643
rlm@2 644 (deftype Symbolic
rlm@2 645 [type
rlm@2 646 expression])
rlm@2 647
rlm@2 648 (defn print-expression [s]
rlm@2 649 (print (.expression s)))
rlm@2 650
rlm@2 651 (defn symbolic-number
rlm@2 652 [symbol]
rlm@2 653 (Symbolic. java.lang.Number symbol))
rlm@2 654
rlm@2 655 (defn simplify-expression [s]
rlm@2 656 (let [e (.expression s)]
rlm@2 657 (cond
rlm@2 658 (not(list? e)) e
rlm@2 659 (= (first e) '+)
rlm@2 660 )
rlm@2 661
rlm@2 662
rlm@2 663
rlm@2 664 (defmethod + [Symbolic Symbolic]
rlm@2 665 [x y]
rlm@2 666 (Symbolic. (.type x) (list '+ (.expression x) (.expression y))))
rlm@2 667
rlm@2 668 (defmethod + [Symbolic java.lang.Number]
rlm@2 669 [s x]
rlm@2 670 (if (zero? x)
rlm@2 671 s
rlm@2 672 (Symbolic. (.type s) (list '+ (.expression s) x))))
rlm@2 673
rlm@2 674 (defmethod sin Symbolic
rlm@2 675 [s]
rlm@2 676 (Symbolic. (.type s) (list 'sin (.expression s))))
rlm@2 677
rlm@2 678 #+end_src
rlm@2 679
rlm@2 680 * COMMENT To-be-tangled Source
rlm@2 681 #+begin_src clojure :tangle utils.clj
rlm@2 682 (ns sicm.utils)
rlm@2 683
rlm@2 684 <<tuples>>
rlm@2 685 <<matrices>>
rlm@2 686 <<arith-tuple>>
rlm@2 687
rlm@2 688 <<differential>>
rlm@2 689 #+end_src
rlm@2 690