Testing Tempo via Direct & Numerical Integration

diffeq.clj

(ns explore_overtone.diffeq)

;; we start with a differential equation:
;;
;; tempo = f(beat) or dx/dt = bps-fn(x)
;;
;; Solving the diff eq via separable eq. gives: (thanks to Tony Grabowski for his help)
;;
;; dt = dx/bps-fn(x) so t = Integrate 1/bps-fn(x) dt or Integrate spb-fn(x) dt
;;  t = log(b + mx) / m - log(b)  (since we set t = 0 at x = 0)
;;
;; Then, algebra gives
;;  x = (b^m)(e^(mt) - b)/m
;;
;; But, we should also still be able to solve the same equation via numerical means:
;;   beat = integrate bps-fn() dt
;; and
;;   time = integrate spb-fn() dBeats

(defn lin
  "dx/dt = m x + b"
  [m b x]
  (+ b (* m x)))

(defn lin-beat2time
  "WolframAlpha: 't = Integrate 1/(mx+b) dx' t = log(b + mx) / m + C.  C=-log(b)"
  [m b x]
  (- (/ (Math/log (+ b (* m x))) m) (Math/log b)))

(defn lin-time2beat
  "given t = log(b + mx) / m - log(b).  Algebra to solve for x = (b^m)(e^(mt) - b)/m.  On WolframAlpha 'Solve ( t = log(b + m*x)/m - log(b)) for x'"
  [m b t]
  (/ (- (* (Math/pow b m) (Math/exp (* m t))) b) m))

;; Integration functions from John Lawrence Aspden blog discussion
;; See: http://www.learningclojure.com/2011/05/numerical-integration-better_26.html
(defn booles-rule [f a b]
  (let [midpoint1 (/ (+ a a a b) 4)
        midpoint2 (/ (+ a a b b) 4)
        midpoint3 (/ (+ a b b b) 4)]
    (* 1/90
       (- b a)
       (+ (* 7 (f a))
          (* 32 (f midpoint1))
          (* 12 (f midpoint2))
          (* 32 (f midpoint3))
          (* 7 (f b))))))

(defn adaptive-rule-recurse [rule f a b desired-error]
  (let [guess (rule f a b)
        midpoint (/ (+ a b) 2)
        better-guess (+ (rule f a midpoint) (rule f midpoint b))
        error-estimate (- guess better-guess)
        abs-error-estimate (if (> error-estimate 0) error-estimate (- error-estimate))]
    (if (< abs-error-estimate desired-error) better-guess
        (let [half-desired-error (/ desired-error 2)]
          (+ (adaptive-rule-recurse rule f a midpoint half-desired-error)
             (adaptive-rule-recurse rule f midpoint b half-desired-error))))))

(defn integrate [f a b]
  (adaptive-rule-recurse booles-rule f a b 0.01))

;; Euler's method http://www.learningclojure.com/2010/09/clojure-is-fast.html
(defn euler [f t0 y0 h its]
  (if (> its 0) 
    (let [t1 (+ t0 h)
          y1 (+ y0 (* h (f y0)))]
      (recur f t1 y1 h (dec its)))
    y0))

(defn do-euler [f t]
  (euler f 0.0 0.0 (/ 1.0 10000) (* 10000 t)))

;; Simple Test #1
;; bps (y) 1.0 -> 2.0 over beat (x) from 0.0 to 100.0
;; bps = (1/100) * x + 1.0 = 0.01 * x + 1.0
(def bps-fn (partial lin 0.01 1.0))
(defn spb-fn [b] (/ 1.0 (bps-fn b)))
(def beat2time1 (partial lin-beat2time 0.01 1.0))
(def time2beat1 (partial lin-time2beat 0.01 1.0))
;; numerical analysis needs two different methods...
(def beat2time2 (partial integrate spb-fn 0.0))
(def time2beat2 (partial do-euler bps-fn))

;; look at it over beats 0, 10, ... 90
(dotimes [x 10]
  (let [cur-beat (* 10.0 x)
        cur-bps (bps-fn cur-beat)
        cur-spb (spb-fn cur-beat)
        cur-time1 (beat2time1 cur-beat)   ;; via direct calc
        cur-beat1 (time2beat1 cur-time1)  ;; via direct calc
        cur-time2 (float (beat2time2 cur-beat))   ;; via numerical integration
        cur-beat2 (float (time2beat2 cur-time2))] ;; via numerical integration
    (println
     (format "%6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f"
             cur-beat cur-bps cur-spb cur-time1 cur-beat1 cur-time2 cur-beat2))
    ))

;; Now, with the proper algorithm (Euler's method), we get
;; expected results from all calculations of time or beat.

;;cur-beat  bps     spb  time1  beat1  time2  beat2
;;  0.000  1.000  1.000  0.000  0.000  0.000  0.000
;; 10.000  1.100  0.909  9.531 10.000  9.531 10.000
;; 20.000  1.200  0.833 18.232 20.000 18.232 20.000
;; 30.000  1.300  0.769 26.236 30.000 26.236 30.000
;; 40.000  1.400  0.714 33.647 40.000 33.647 40.000
;; 50.000  1.500  0.667 40.547 50.000 40.547 50.000
;; 60.000  1.600  0.625 47.000 60.000 47.000 60.000
;; 70.000  1.700  0.588 53.063 70.000 53.063 70.000
;; 80.000  1.800  0.556 58.779 80.000 58.779 80.000
;; 90.000  1.900  0.526 64.185 90.000 64.185 90.000