Executable examples of the SICM book

(note on first visit: this page is slow to load and maybe needs reload after some seconds)

This html-file covers the first chapter of the SICM book [1], using the great Emmy computer algebra system [2]. The examples are also available via visual programming [3]. Further details can be found in part 7 of my blog series [4].

[1] SICM book (unofficial version with nice rendering)

[2] Emmy

[3] Visual programming with clj-tiles

[4] Blog post

This html-file is hand-written (i.e. not generated from a Jupyter notebook or the like). All code on the page can be changed to immediate effect, thanks to Klipse's browser execution. The SICM examples are calculated using a tailor-made Scheme syntax to JavaScript compiler of only about 160 lines of code.

This html file, being the sole source code for both the examples and their syntax compiler, is meant to also be modified and run locally. Just copy-paste the html-code into your text-editor (e.g. Notepad++), save locally and open in your browser.

Some Code Snippets


1 + 2 + 3;
    

#t
    

The Compiler


var symbol = emmy.symbol;

var specialchars =
[["", /\:/g],
 ["__gt_", /->/g],
 ["symbol_$1", /\'(\w+)/g],
 ["$1_dot_$2", /(\d+)\.(\d+)/g],
 [" minus_$2", /(\s+)\-(\w+)/g],
 ["$1_", /(\w+)\-/g],
 ["[ / $1 $2 ],", /(\w+)\/(\d+)/g]];

var post_specialchars =
[["$1.$2", /(\d+)_dot_(\d+)/g],
 ["$1e-$2", /(\d+)e_(\d+)/],
 ["-$1", /minus_(\w+)/g],
 //case symbol with hat 'v_r^x
 ["symbol_$1^$2", /symbol_(\w+)expt(\w+)/]]

var mathfns =
[["div", /\//g],
["mul", /\*/g],
["sub", /\-/g],
["add", /\+/g],
["expt", /\^/g]]

mathfns.map(x=> x[0]).map(loadEnv);

var replaceMath = (txt) =>
specialchars.concat(mathfns)
.reduce((s,r) => s.replace(r[1],r[0]), txt);

var swapFirst = (j) =>
j.constructor == Array && mathfns.map(x=> x[0]).includes(j[1])
? [j[1], j[0]].concat(j.slice(2)).map(swapFirst)
: j.constructor == Array
? j.map(swapFirst)
: j;

var insertCommas = (txt) =>
txt.replace(/(\w+)/g,'"$1",')
.replace(/\,\s+\]/g," ]");

var makeBrackets = (txt) =>
txt.trim()
.replace(/\(/g,"[ ")
.replace(/\)/g," ],")
.replace(/,$/,"");

var textToJson = (txt) =>
insertCommas(replaceMath(makeBrackets(txt)));

var drv = emmy.D;

var replaceD = (j) =>
j.constructor == Array
? j.map(replaceD)
: j.constructor == String && j == "D"
? "drv"
: j.constructor == String
? post_specialchars.reduce((s,r) => s.replace(r[1],r[0]), j)
: j;

["to_infix", "simplify"].map(loadEnv);

var preAmble = (j) =>
j[0] != "define"
? ["to_infix", j]
: j;

var identity = x => x;

var edgeCases = (j) =>
j[0] === "define"
&& j[2].constructor == Array && j[2][0] === "let"
&& j[1].constructor == String
? ["identity", "'let not allowed in variable definition'"]
: j;

var letFlat2 = (j) =>
j[0] === "let"
? j[1].concat(letFlat2(j[2]))
: [j];

var constructLet = (j) =>
["let", j.slice(0, j.length - 1), j[j.length - 1]]

var letFlat = (j) =>
j.constructor == Array && j[0] === "let"
? constructLet(letFlat2(j))
: j.constructor == Array
? j.map(letFlat)
: j;

var modifyJson = (j) =>
preAmble(replaceD(letFlat(edgeCases(swapFirst(j)))));

var makeFun = (j, callBack) =>
"var " + j[1][0] + " = (" + j[1].slice(1) + ") => "
+ callBack(j[2]) +";";

var makeFunFun = (j, callBack) =>
"var " + j[1][0][0] + " = (" + j[1][0].slice(1) + ") => "
+ "(" + j[1].slice(1) +") => " + callBack(j[2]) + ";";

var smap = (f, v) =>
v.length === 1
? f(v[0]) +";"
: f(v[0]) +"; "  + smap(f, v.slice(1));

var makeLet = (j, callBack) =>
"{" + smap((l) => "let " + l[0] + " = " + callBack(l[1]), j[1])
+ " return " +  callBack(j[2]) + "; }"

var jsonToJs = (j) =>
j.constructor == Array && j[0] === "let"
? makeLet(j, jsonToJs)
:j.constructor == Array && j[0] === "define"
&& j[1].constructor == Array
&& j[1][0].constructor == Array
? makeFunFun(j, jsonToJs)
:j.constructor == Array && j[0] === "define"
&& j[1].constructor == Array
? makeFun(j, jsonToJs)
: j.constructor == Array && j[0] === "lambda"
? "(" + j[1] + ") => " + jsonToJs(j[2])
:j.constructor == Array && j[0] === "define"
&& j[1].constructor == String
? "var " + j[1] + " = " + jsonToJs(j[2]) + ";"
:j.constructor == String && j.substring(0, 7) == "symbol_"
? 'symbol("' + j.substring(7, j.length) + '")'
:j.constructor == Array
? jsonToJs(j[0]) + "("  +  j.slice(1).map(jsonToJs) + ")"
: j;

var expressionToJs = (expr) =>
jsonToJs(modifyJson(JSON.parse(textToJson(expr))));

["sin", "cos", "pi", "velocity", "dot_product",
"up", "nth", "Gamma", "compose", "literal_function",
"definite_integral", "coordinate", "minimize",
"linear_interpolants", "multidimensional_minimize",
"make_path", "count",
"square", "partial", "ref", "sqrt", "atan", "down",
"solve_linear_left", "state_advancer", "cross_product",
"Rx", "Ry", "Rz", "osculating_path"].map(loadEnv);

var show_expression = simplify;
var vector_length = count;
var velocities = velocity;
var coordinates = coordinate;
var time = state => nth(state, 0);
var get = (obj, field) => obj[field.toString()]; null;
    

Warm up


(+ 1 2 3) 
    

(cos :pi) 
    

(up (* 'R (cos 'phi))
    (* 'R (sin 'phi))) 
    

Remember that the generated Javascript is shown by adding a semicolon at the beginning


;(up (* 'R (cos 'phi))
    (* 'R (sin 'phi)))
    

1 Lagrangian Mechanics

1.4 Computing Actions


(define ((L-free-particle mass) local)
(let ((v (velocity local)))
(* 1/2 mass (dot-product v v)))) 
    

(define q
 (up (literal-function 'x)
 (literal-function 'y)
 (literal-function 'z))) 
    

((Gamma q) 't) 
    

((compose (L-free-particle 'm) (Gamma q)) 't) 
    

(define (Lagrangian-action L q t1 t2)
  (definite-integral (compose L (Gamma q)) t1 t2)) 
    

(define (test-path t)
  (up (+ (* 4 t) 7)
      (+ (* 3 t) 5)
      (+ (* 2 t) 1))) 
    

(Lagrangian-action (L-free-particle 3.0)
                   test-path 0.0 10.0) 
    

(define ((make-eta nu t1 t2) t)
(* (- t t1) (- t t2) (nu t))) 
    

(define ((varied-free-particle-action mass q nu t1 t2) eps)
  (let ((eta (make-eta nu t1 t2)))
    (Lagrangian-action (L-free-particle mass)
                       (+ q (* eps eta))
                       t1
                       t2))) 
    

((varied-free-particle-action 3.0 test-path
                              (up sin cos square)
                              0.0 10.0)
 0.001) 
    

Calculation takes a few seconds. Add a blank at the end to start.


(get
 (minimize
  (varied-free-particle-action 3.0 test-path
   (up sin cos square)
    0.0 10.0)
  -2.0 1.0) 'value)
    

((make-path 0 0 10 5 (up -1 2 -3 4)) 6) 
    

(define ((parametric-path-action Lagrangian t0 q0 t1 q1) qs)
  (let ((path (make-path t0 q0 t1 q1 qs)))
    (Lagrangian-action Lagrangian path t0 t1))) 
    

(define (find-path Lagrangian t0 q0 t1 q1 n)
  (let ((initial-qs (linear-interpolants q0 q1 n)))
    (let ((minimizing-qs
            (multidimensional-minimize
              (parametric-path-action Lagrangian t0 q0 t1 q1)
              initial-qs)))
      (make-path t0 q0 t1 q1 minimizing-qs)))) 
    

(define ((L-harmonic m k) local)
  (let ((q (coordinate local))
        (v (velocity local)))
    (- (* 1/2 m (square v)) (* 1/2 k (square q))))) 
    

Calculation takes a few seconds, add a blank at the and to start.


(define q-harmonic 
  (find-path (L-harmonic 1.0 1.0) 0.0 1.0 :pi/2 0.0 3))
    

Add a blank at the and to start


(- (cos 0.8) (q-harmonic 0.8))
    

1.5 The Euler–Lagrange Equations

1.5.2 Computing Lagrange's Equations


(define ((Lagrange-equations Lagrangian) q)
  (- (D (compose ((partial 2) Lagrangian) (Gamma q)))
     (compose ((partial 1) Lagrangian) (Gamma q)))) 
    

Here I had to rename test-path to general-test-path to avoid, in chapter 1.4, conflict with the minimize procedure there.


(define (general-test-path t)
  (up (+ (* 'a t) 'a0)
      (+ (* 'b t) 'b0)
      (+ (* 'c t) 'c0))) 
    

(((Lagrange-equations (L-free-particle 'm))
  general-test-path)
 't) 
    

(show-expression
  (((Lagrange-equations (L-free-particle 'm))
    (literal-function 'x))
   't)) 
    

(define (proposed-solution t)
  (* 'A (cos (+ (* 'omega t) 'phi)))) 
    

(show-expression
  (((Lagrange-equations (L-harmonic 'm 'k))
    proposed-solution)
   't)) 
    

Exercise 1.11: Kepler's third law

Show that a planet in circular orbit satisfies Kepler's third law n^2a^3=G(M_1 + m_2), where n is the angular frequency of the orbit and a is the distance between sun and planet. (Hint: use the reduced mass to construct the Lagrangian)


(define ((L-Kepler-central-polar m V) local)
  (let ((q (coordinate local))
        (qdot (velocity local)))
    (let ((r (ref q 0))     (phi (ref q 1))
                            (rdot (ref qdot 0)) (phidot (ref qdot 1)))
      (- (* 1/2 m
            (+ (square rdot) (square (* r phidot))) )
         (V r))))) 
    

(define ((gravitational-energy G m1 m2) r)
  (- (/ (* G m1 m2) r))) 
    

(define (circle t)
  (up 'a (* 'n t))) 
    

(define lagrangian-reduced
(L-Kepler-central-polar (/ (* 'M_1 'm_2) (+ 'M_1 'm_2))
(gravitational-energy 'G 'M_1 'm_2))) 
    

(show-expression
(((Lagrange-equations lagrangian-reduced) circle) 't)) 
    

1.6 How to find Lagrangians


(define ((L-uniform-acceleration m g) local)
  (let ((q (coordinate local))
        (v (velocity local)))
    (let ((y (ref q 1)))
      (- (* 1/2 m (square v)) (* m g y))))) 
    

(show-expression
  (((Lagrange-equations
      (L-uniform-acceleration 'm 'g))
    (up (literal-function 'x)
        (literal-function 'y)))
   't)) 
    

(define ((L-central-rectangular m U) local)
  (let ((q (coordinate local))
        (v (velocity local)))
    (- (* 1/2 m (square v))
       (U (sqrt (square q)))))) 
    

(show-expression
  (((Lagrange-equations
      (L-central-rectangular 'm (literal-function 'U)))
    (up (literal-function 'x)
        (literal-function 'y)))
   't)) 
    

(show-expression
  (((Lagrange-equations
      (L-Kepler-central-polar 'm (literal-function 'U)))
    (up (literal-function 'r)
        (literal-function 'phi)))
   't)) 
    

1.6.1 Coordinate Transformations


(define ((F->C F) local)
  (up (time local)
      (F local)
      (+ (((partial 0) F) local)
         (* (((partial 1) F) local)
            (velocity local))))) 
    

(define (p->r local)
  (let ((polar-tuple (coordinate local)))
    (let ((r (ref polar-tuple 0))
          (phi (ref polar-tuple 1)))
      (let ((x (* r (cos phi)))
            (y (* r (sin phi))))
        (up x y))))) 
    

(show-expression
  (velocity
    ((F->C p->r)
     (up 't (up 'r 'phi) (up 'rdot 'phidot))))) 
    

(define (L-central-polar m U)
  (compose (L-central-rectangular m U) (F->C p->r))) 
    

(show-expression
  ((L-central-polar 'm (literal-function 'U))
   (up 't (up 'r 'phi) (up 'rdot 'phidot)))) 
    

Coriolis and centrifugal forces


(define ((L-free-rectangular m) local)
  (let ((vx (ref (velocities local) 0))
        (vy (ref (velocities local) 1)))
    (* 1/2 m (+ (square vx) (square vy))))) 
    

(define (L-free-polar m)
  (compose (L-free-rectangular m) (F->C p->r))) 
    

(define ((F Omega) local)
  (let ((t (time local))
        (r (ref (coordinates local) 0))
        (theta (ref (coordinates local) 1)))
    (up r (+ theta (* Omega t))))) 
    

(define (L-rotating-polar m Omega)
  (compose (L-free-polar m) (F->C (F Omega)))) 
    

(define (L-rotating-rectangular m Omega)
  (compose (L-rotating-polar m Omega) (F->C r->p))) 
    

r->p added


(define (r->p local)
  (let ((rect-tuple (coordinate local)))
    (let ((x (ref rect-tuple 0))
          (y (ref rect-tuple 1)))
       (let ((r (sqrt (square rect-tuple)))
             (phi (atan (/ y x))))
         (up r phi))))) 
    

(show-expression
((L-rotating-rectangular 'm 'Omega)
(up 't (up 'x_r 'y_r) (up 'xdot_r 'ydot_r)))) 
    

(+ (* 1/2 (expt 'Omega 2) 'm (expt 'x_r 2))
(* 1/2 (expt 'Omega 2) 'm (expt 'y_r 2))
(* -1 'Omega 'm 'xdot_r 'y_r)
(* 'Omega 'm 'ydot_r 'x_r)
(* 1/2 'm (expt 'xdot_r 2))
(* 1/2 'm (expt 'ydot_r 2))) 
    

x_r, y_r: underscore added. Calculation takes a few seconds, add a blank at the and to start


(show-expression
(((Lagrange-equations (L-rotating-rectangular 'm 'Omega))
  (up (literal-function 'x_r) (literal-function 'y_r)))
 't))
    

definitions x_r y_r added


(define x_r (literal-function 'x_r)) 
    

(define y_r (literal-function 'y_r)) 
    

(down
(+ (* -1 (expt 'Omega 2) 'm (x_r 't))
(* -2 'Omega 'm ((D y_r) 't))
(* 'm (((expt D 2) x_r) 't)))
(+ (* -1 (expt 'Omega 2) 'm (y_r 't))
(* 2 'Omega 'm ((D x_r) 't))
(* 'm (((expt D 2) y_r) 't)))) 
    

1.6.2 Systems with Rigid Constraints

A pendulum driven at the pivot

See here for a presentation of the Driven Pendulum using visual programming


(define ((T-pend m l g ys) local)
(let ((t (time local))
(theta (coordinate local))
(thetadot (velocity local)))
(let ((vys (D ys)))
(* 1/2 m
(+ (square (* l thetadot))
(square (vys t))
(* 2 l (vys t) thetadot (sin theta))))))) 
    

(define ((V-pend m l g ys) local)
(let ((t (time local))
      (theta (coordinate local)))
  (* m g (- (ys t) (* l (cos theta)))))) 
    

Because used later, rename L-pend to L-pendulum


(define L-pendulum (- T-pend V-pend)) 
    

(show-expression
(((Lagrange-equations
(L-pendulum 'm 'l 'g (literal-function 'y_s)))
(literal-function 'theta))
't)) 
    

1.6.3 Constraints as Coordinate Transformations

L-uniform-acceleration already defined above


(define ((L-uniform-acceleration m g) local)
  (let ((q (coordinate local))
        (v (velocity local)))
     (let ((y (ref q 1)))
       (- (* 1/2 m (square v)) (* m g y))))) ;;execution not necessary
    

(define ((dp-coordinates l y_s) local)
(let ((t (time local))
(theta (coordinate local)))
(let ((x (* l (sin theta)))
(y (- (y_s t) (* l (cos theta)))))
(up x y)))) 
    

(define (L-pend m l g y_s)
(compose (L-uniform-acceleration m g)
(F->C (dp-coordinates l y_s)))) 
    

(show-expression
((L-pend 'm 'l 'g (literal-function 'y_s))
(up 't 'theta 'thetadot))) 
    

1.7 Evolution of Dynamical State


(define (Lagrangian->acceleration L)
(let ((P ((partial 2) L)) (F ((partial 1) L)))
(solve-linear-left
((partial 2) P)
(- F
(+ ((partial 0) P)
(* ((partial 1) P) velocity)))))) 
    

(define (Lagrangian->state-derivative L)
(let ((acceleration (Lagrangian->acceleration L)))
(lambda (state)
(up 1
(velocity state)
(acceleration state))))) 
    

(define (harmonic-state-derivative m k)
(Lagrangian->state-derivative (L-harmonic m k))) 
    

(show-expression
((harmonic-state-derivative 'm 'k)
(up 't (up 'x 'y) (up 'v_x 'v_y)))) 
    

(up 1 (up 'v_x 'v_y) (up (/ (* -1 'k 'x) 'm) (/ (* -1 'k 'y) 'm))) 
    

(define ((Lagrange-equations-first-order L) q v)
  (let ((state-path (qv->state-path q v)))
    (- (D state-path)
       (compose (Lagrangian->state-derivative L)
                state-path)))) 
    

(define ((qv->state-path q v) t)
  (up t (q t) (v t))) 
    

(show-expression
 (((Lagrange-equations-first-order (L-harmonic 'm 'k))
   (up (literal-function 'x)
       (literal-function 'y))
   (up (literal-function 'v_x)
       (literal-function 'v_y)))
  't)) 
    

Numerical integration


((state-advancer harmonic-state-derivative 2.0 1.0)
(up 1.0 (up 1.0 2.0) (up 3.0 4.0))
10.0
1.0e-12)
    

(up 11.0
    (up 3.7127916645844437 5.420620823651583)
    (up 1.6148030925459782 1.8189103724750855)) 
    

(define ((periodic-drive amplitude frequency phase) t)
(* amplitude (cos (+ (* frequency t) phase)))) 
    

(define (L-periodically-driven-pendulum m l g A omega)
(let ((ys (periodic-drive A omega 0)))
(L-pend m l g ys))) 
    

(show-expression
(((Lagrange-equations
(L-periodically-driven-pendulum 'm 'l 'g 'A 'omega))
(literal-function 'theta))
't)) 
    

(define (pend-state-derivative m l g A omega)
(Lagrangian->state-derivative
(L-periodically-driven-pendulum m l g A omega))) 
    

(show-expression
((pend-state-derivative 'm 'l 'g 'A 'omega)
(up 't 'theta 'thetadot))) 
    

1.8 Conserved Quantities

1.8.2 Energy Conservation


(define (Lagrangian->energy L)
(let ((P ((partial 2) L)))
(- (* P velocity) L))) 
    

1.8.3 Central Forces in Three Dimensions


(define ((T3-spherical m) state)
(let ((q (coordinate state))
(qdot (velocity state)))
(let ((r (ref q 0))
(theta (ref q 1))
(rdot (ref qdot 0))
(thetadot (ref qdot 1))
(phidot (ref qdot 2)))
(* 1/2 m
(+ (square rdot)
(square (* r thetadot))
(square (* r (sin theta) phidot))))))) 
    

Change the second define into a let


(define (L3-central m Vr)
  (let ((Vs (lambda (state)
               (let ((r (ref (coordinate state) 0)))
                 (Vr r)))))
     (- (T3-spherical m) Vs))) 
    

(show-expression
(((partial 1) (L3-central 'm (literal-function 'V)))
(up 't
(up 'r 'theta 'phi)
(up 'rdot 'thetadot 'phidot)))) 
    

(show-expression
(((partial 2) (L3-central 'm (literal-function 'V)))
(up 't
(up 'r 'theta 'phi)
(up 'rdot 'thetadot 'phidot)))) 
    

(define ((ang-mom-z m) rectangular-state)
(let ((xyz (coordinate rectangular-state))
(v (velocity rectangular-state)))
(ref (cross-product xyz (* m v)) 2))) 
    

(define (s->r spherical-state)
(let ((q (coordinate spherical-state)))
(let ((r (ref q 0))
(theta (ref q 1))
(phi (ref q 2)))
(let ((x (* r (sin theta) (cos phi)))
(y (* r (sin theta) (sin phi)))
(z (* r (cos theta))))
(up x y z))))) 
    

(show-expression
((compose (ang-mom-z 'm) (F->C s->r))
(up 't
(up 'r 'theta 'phi)
(up 'rdot 'thetadot 'phidot)))) 
    

(show-expression
((Lagrangian->energy (L3-central 'm (literal-function 'V)))
(up 't
(up 'r 'theta 'phi)
(up 'rdot 'thetadot 'phidot)))) 
    

1.8.4 The Restricted Three-Body Problem


(define ((L0 m V) local)
(let ((t (time local))
(q (coordinates local))
(v (velocities local)))
(- (* 1/2 m (square v)) (V t q)))) 
    

(define ((V a GM0 GM1 m) t xy)
(let ((Omega (sqrt (/ (+ GM0 GM1) (expt a 3))))
(a0 (* (/ GM1 (+ GM0 GM1)) a))
(a1 (* (/ GM0 (+ GM0 GM1)) a)))
(let ((x (ref xy 0)) (y (ref xy 1))
(x0 (* -1 a0 (cos (* Omega t))))
(y0 (* -1 a0 (sin (* Omega t))))
(x1 (* +1 a1 (cos (* Omega t))))
(y1 (* +1 a1 (sin (* Omega t)))))
(let ((r0
(sqrt (+ (square (- x x0)) (square (- y y0)))))
(r1
(sqrt (+ (square (- x x1)) (square (- y y1))))))
(- (+ (/ (* GM0 m) r0) (/ (* GM1 m) r1))))))) 
    

(define ((LR3B m a GM0 GM1) local)
(let ((q (coordinates local))
(qdot (velocities local))
(Omega (sqrt (/ (+ GM0 GM1) (expt a 3))))
(a0 (* (/ GM1 (+ GM0 GM1)) a))
(a1 (* (/ GM0 (+ GM0 GM1)) a)))
(let ((x (ref q 0))     (y (ref q 1))
(xdot (ref qdot 0)) (ydot (ref qdot 1)))
(let ((r0 (sqrt (+ (square (+ x a0)) (square y))))
(r1 (sqrt (+ (square (- x a1)) (square y)))))
(+ (* 1/2 m (square qdot))
(* 1/2 m (square Omega) (square q))
(* m Omega (- (* x ydot) (* xdot y)))
(/ (* GM0 m) r0) (/ (* GM1 m) r1)))))) 
    

(define ((LR3B1 m a0 a1 Omega GM0 GM1) local)
(let ((q (coordinates local))
(qdot (velocities local)))
(let ((x (ref q 0))     (y (ref q 1))
(xdot (ref qdot 0)) (ydot (ref qdot 1)))
(let ((r0 (sqrt (+ (square (+ x a0)) (square y))))
(r1 (sqrt (+ (square (- x a1)) (square y)))))
(+ (* 1/2 m (square qdot))
(* 1/2 m (square Omega) (square q))
(* m Omega (- (* x ydot) (* xdot y)))
(/ (* GM0 m) r0) (/ (* GM1 m) r1)))))) 
    

((Lagrangian->energy (LR3B1 'm 'a_0 'a_1 'Omega 'GM_0 'GM_1))
(up 't (up 'x_r 'y_r) (up 'v_r^x 'v_r^y))) 
    

(+ (* 1/2 'm (expt 'v_r^x 2))
(* 1/2 'm (expt 'v_r^y 2))
(/ (* -1 'GM_0 'm)
(sqrt (+ (expt (+ 'x_r 'a_0) 2) (expt 'y_r 2))))
(/ (* -1 'GM_1 'm)
(sqrt (+ (expt (- 'x_r 'a_1) 2) (expt 'y_r 2))))
(* -1/2 'm (expt 'Omega 2) (expt 'x_r 2))
(* -1/2 'm (expt 'Omega 2) (expt 'y_r 2))) 
    

1.8.5 Noether’s Theorem


(define (F-tilde angle-x angle-y angle-z)
(compose (Rx angle-x) (Ry angle-y) (Rz angle-z) coordinate)) 
    

L-central-rectangular is already defined in chapter 1.6


(define ((L-central-rectangular m U) state)
  (let ((q (coordinate state))
        (v (velocity state)))
    (- (* 1/2 m (square v))
       (U (sqrt (square q)))))) ;;execution not needed
    

A let within a variable definition is not allowed in our little Scheme compiler,


(define the-Noether-integral
  (let ((L (L-central-rectangular
              'm (literal-function 'U))))
     (* ((partial 2) L) ((D F-tilde) 0 0 0)))) 
    

... so we split in two expressions.


(define let-L (L-central-rectangular
              'm (literal-function 'U))) 
    

(define the-Noether-integral
  (* ((partial 2) let-L) ((D F-tilde) 0 0 0))) 
    

(show-expression
(the-Noether-integral
(up 't
(up 'x 'y 'z)
(up 'vx 'vy 'vz)))) 
    

(down (+ (* 'm 'vy 'z) (* -1 'm 'vz 'y))
(+ (* 'm 'vz 'x) (* -1 'm 'vx 'z))
(+ (* 'm 'vx 'y) (* -1 'm 'vy 'x))) 
    

1.9 Abstraction of Path Functions


(define ((Gamma-bar f-bar) local)
((f-bar (osculating-path local)) (time local))) 
    

define withing define is not allowed in our little Scheme compiler ...


(define (F->C1 F)
(define (C local)
(let ((n (vector-length local)))
(define (f-bar q-prime)
(define q
(compose F (Gamma q-prime)))
(Gamma q n))
((Gamma-bar f-bar) local)))
C) ;;do not execute
    

... so we use let. Moreover, F->C is already defined in chapter 1.6.1, so we rename to F->C1


(define (F->C1 F)
  (let ((C (lambda (local)
    (let ((n (vector-length local)  )
          (f-bar (lambda (q-prime)
                   (let ((q (compose F (Gamma q-prime))  ))
                     (Gamma q n)))))
      ((Gamma-bar f-bar) local)))))
    C)) 
    

(show-expression
((F->C1 p->r)
(up 't (up 'r 'theta) (up 'rdot 'thetadot)))) 
    

define withing define is not allowed in our little Scheme compiler ...


(define (Dt F)
(define (DtF state)
(let ((n (vector-length state)))
(define (DF-on-path q)
(D (compose F (Gamma q (- n 1)))))
((Gamma-bar DF-on-path) state)))
DtF) ;;do not execute
    

... so we use let


(define (Dt F)
  (let ((DtF (lambda (state)
               (let ((n (vector-length state))
                   (DF-on-path (lambda (q)
                                 (D (compose F (Gamma q (- n 1)))))))
                 ((Gamma-bar DF-on-path) state)))))
    DtF)) 
    

(define (Euler-Lagrange-operator L)
(- (Dt ((partial 2) L)) ((partial 1) L))) 
    

((Euler-Lagrange-operator
   (L-harmonic 'm 'k))
     (up 't 'x 'v 'a)) 
    

(+ (* 'a 'm) (* 'k 'x)) 
    

((compose
(Euler-Lagrange-operator (L-harmonic 'm 'k))
(Gamma (literal-function 'x) 4))
't) 
    

(+ (* 'k ((literal-function 'x) 't))
   (* 'm (((expt D 2) (literal-function 'x)) 't)))