sicmutils · sritchie · Jan 13, 2022 · Jan 13, 2022 · Jan 13, 2022 · Jan 13, 2022
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -2,6 +2,21 @@
 
 ## [unreleased]
 
+- #456:
+
+  - `sicmutils.mechanics.lagrange/{Γ,Γ-bar}` are removed in favor of the
+    existing `Gamma` and `Gamma-bar` functions. The `sicmutils.env` aliases are
+    gone as well.
+
+  - `sicmutils.mechanics.lagrange/Lagrange-interpolation-function` now returns
+    an actual polynomial instance. Because polynomials support `IFn` and respond
+    to the derivative operator `D`, this makes the `find-path` example on pages
+    22/23 of SICM run about 5x faster.
+
+  - Richardson extrapolation is now implemented as a functional fold. The
+    exposition in `sicmutils.polynomial.richardson` discusses this; the
+    namespaces gains `richardson-fold`, `richardson-sum` and `richardson-scan`.
+
 - #455 makes `sicmutils.util.aggregate/scan` and
   `sicmutils.algebra.fold/fold->scan-fn` slightly more efficient by dropping the
   first element of the returned sequence before mapping the `present` function.

diff --git a/src/sicmutils/env.cljc b/src/sicmutils/env.cljc
@@ -461,8 +461,8 @@ constant [Pi](https://en.wikipedia.org/wiki/Pi)."}
   ->local
   Euler-Lagrange-operator
   F->C
-  Gamma Γ
-  Gamma-bar Γ-bar
+  Gamma
+  Gamma-bar
   Lagrange-equations
   Lagrange-equations-first-order
   Lagrange-interpolation-function

diff --git a/src/sicmutils/mechanics/lagrange.cljc b/src/sicmutils/mechanics/lagrange.cljc
@@ -24,6 +24,7 @@
             [sicmutils.calculus.derivative :refer [D partial]]
             [sicmutils.function :as f]
             [sicmutils.generic :as g :refer [cos sin + - * /]]
+            [sicmutils.polynomial :as p]
             [sicmutils.structure :refer [up down]]
             [sicmutils.function :as f :refer [compose]]))
 
@@ -137,47 +138,34 @@
      (fn [t]
        (->> Dqs (map #(% t)) (cons t) (apply up))))))
 
-(def Γ Gamma)
+(defn Lagrangian-action [L q t1 t2]
+  (q/definite-integral
+    (compose L (Gamma q)) t1 t2
+    {:compile? false}))
 
-(defn Lagrangian-action
-  [L q t1 t2]
-  (q/definite-integral (compose L (Γ q)) t1 t2 {:compile? false}))
-
-(defn Lagrange-equations
-  [Lagrangian]
+(defn Lagrange-equations [Lagrangian]
   (fn [q]
-    (- (D (compose ((partial 2) Lagrangian) (Γ q)))
-       (compose ((partial 1) Lagrangian) (Γ q)))))
+    (- (D (compose ((partial 2) Lagrangian) (Gamma q)))
+       (compose ((partial 1) Lagrangian) (Gamma q)))))
 
-(defn linear-interpolants
-  [x0 x1 n]
+(defn linear-interpolants [x0 x1 n]
   (let [n+1 (inc n)
-        dx (/ (- x1 x0) n+1)]
+        dx  (/ (- x1 x0) n+1)]
     (for [i (range 1 n+1)]
       (+ x0 (* i dx)))))
 
 (defn Lagrange-interpolation-function
+  "Given `ys` (a sequence of function values) and `xs` (an equal-length sequence
+  of function inputs), returns a [[sicmutils.polynomial/Polynomial]] instance
+  guaranteed to pass through all supplied `xs` and `ys`.
+
+  The contract for inputs is that `(map vector xs ys)` should return a sequence
+  of pairs of points."
   [ys xs]
-  (let [n (count ys)]
-    (assert (= (count xs) n))
-    (-> (fn [x]
-          (reduce + 0
-                  (for [i (range n)]
-                    (/ (reduce * 1
-                               (for [j (range n)]
-                                 (if (= j i)
-                                   (nth ys i)
-                                   (- x (nth xs j)))))
-                       (let [xi (nth xs i)]
-                         (reduce * 1
-                                 (for [j (range n)]
-                                   (cond (< j i) (- (nth xs j) xi)
-                                         (= j i) (if (odd? i) -1 1)
-                                         :else (- xi (nth xs j))))))))))
-        (f/with-arity [:exactly 1]))))
-
-(defn Lagrangian->acceleration
-  [L]
+  (p/from-points
+   (map vector xs ys)))
+
+(defn Lagrangian->acceleration [L]
   (let [P ((partial 2) L)
         F ((partial 1) L)]
     (/ (- F
@@ -197,16 +185,14 @@
   [q v]
   #(up % (q %) (v %)))
 
-(defn Lagrange-equations-first-order
-  [L]
+(defn Lagrange-equations-first-order [L]
   (fn [q v]
     (let [state-path (qv->state-path q v)]
       (- (D state-path)
          (compose (Lagrangian->state-derivative L)
                   state-path)))))
 
-(defn Lagrangian->energy
-  [L]
+(defn Lagrangian->energy [L]
   (let [P ((partial 2) L)]
     (- (* P velocity) L)))
 
@@ -227,13 +213,10 @@
                    (+ sum (* (nth state0 n) dt-n:n!))
                    (/ (* dt-n:n! dt) n))))))))
 
-(defn Gamma-bar
-  [f]
+(defn Gamma-bar [f]
   (fn [local]
     ((f (osculating-path local)) (first local))))
 
-(def Γ-bar Gamma-bar)
-
 (defn F->C
   "Accepts a coordinate transformation `F` from a local tuple to a new coordinate
   structure, and returns a function from `local -> local` that applies the
@@ -249,13 +232,13 @@
       ((Gamma-bar f-bar) local))))
 
 (defn Dt [F]
-  (let [G-bar (fn [q]
-                (D (compose F (Γ q))))]
-    (Γ-bar G-bar)))
+  (letfn [(G-bar [q]
+            (D (compose F (Gamma q))))]
+    (Gamma-bar G-bar)))
 
-(defn Euler-Lagrange-operator
-  [L]
-  (- (Dt ((partial 2) L)) ((partial 1) L)))
+(defn Euler-Lagrange-operator [L]
+  (- (Dt ((partial 2) L))
+     ((partial 1) L)))
 
 (defn L-rectangular
   "Lagrangian for a point mass on with the potential energy V(x, y)"

diff --git a/src/sicmutils/mechanics/rigid.cljc b/src/sicmutils/mechanics/rigid.cljc
@@ -55,11 +55,11 @@
 
 (defn M->omega
   [M-of-q]
-  (-> M-of-q M-of-q->omega-of-t L/Γ-bar))
+  (-> M-of-q M-of-q->omega-of-t L/Gamma-bar))
 
 (defn M->omega-body
   [M-of-q]
-  (-> M-of-q M-of-q->omega-body-of-t L/Γ-bar))
+  (-> M-of-q M-of-q->omega-body-of-t L/Gamma-bar))
 
 (defn Euler-state->omega-body
   [[_ [θ _ ψ] [θdot φdot ψdot]]]

diff --git a/src/sicmutils/polynomial.cljc b/src/sicmutils/polynomial.cljc
@@ -460,6 +460,14 @@
   (map #(identity n %)
        (range 0 n)))
 
+(defn from-points
+  "Given a sequence of points of the form `[x, f(x)]`, returns a univariate
+  polynomial that passes through each input point.
+
+  The degree of the returned polynomial is equal to `(dec (count xs))`."
+  [xs]
+  (pi/lagrange xs (identity)))
+
 (declare add)
 
 (defn linear
@@ -1599,17 +1607,6 @@
          poly     (x/evaluate expr sym->var operator-table)]
      (cont poly sorted))))
 
-(defn from-points
-  "Given a sequence of points of the form `[x, f(x)]`, returns a univariate
-  polynomial that passes through each input point.
-
-  The degree of the returned polynomial is equal to `(dec (count xs))`."
-  [xs]
-  (expression->
-   (g/simplify
-    (pi/lagrange xs 'x))
-   (fn [p _] p)))
-
 (let [* (sym/symbolic-operator '*)
       + (sym/symbolic-operator '+)
       expt (sym/symbolic-operator 'expt)]

diff --git a/src/sicmutils/polynomial/interpolate.cljc b/src/sicmutils/polynomial/interpolate.cljc
@@ -52,7 +52,7 @@
                                     (get-in points [j 0]))
                            p (apply g/* (map #(g/- x %) others))
                            q (apply g/* (map #(g/- a %) others))]
-                       (g// (g/* fa p) q)))]
+                       (g/* (g/invert q) fa p)))]
     (transduce (map-indexed build-term)
                g/+
                points)))
@@ -580,12 +580,17 @@
 ;; This method reverses the order of the points, since rows are built from the
 ;; bottom up:
 ;;
-;; p4 p34 p234 p1234 p01234
-;; p3 p23 p123 p0123 .
-;; p2 p12 p012 .     .
-;; p1 p01 .    .     .
+;; p4 p43 p432 p4321 p43210
+;; p3 p32 p321 p3210 .
+;; p2 p21 p210 .     .
+;; p1 p10 .    .     .
 ;; p0 .   .    .     .
 ;;
+;; The order of the points is reversed, but this is fine; polynomial
+;; interpolation doesn't care about the order of points. (NOTE that this WILL be
+;; something we have to consider in the fold version of Richardson
+;; extrapolation, in `sicmutils.polynomial.richardson`!)
+;;
 ;; Notice that the /diagonal/ of this tableau is identical to the top row of the
 ;; tableau before the points were reversed.
 ;;
@@ -819,6 +824,3 @@
 ;; - `richardson.cljc` for a specialized implementation of polynomial
 ;;   interpolation, when you know something about the ratios between successive
 ;;   `x` elements in the point sequence.
-;;
-;; NOTE: For bonus points, see if you can figure out how to write Richardson
-;; extrapolation as a functional fold!