Executable examples of the SICM book part 3

This html-file is part of a series of pages [1]. It covers the third part of the SICM book [2].

[1] Html series

[2] SICM book

1 + 2 + 3;
    

#t
    

The Compiler

["pi", "_pi", "state__GT_t", "symbol"].map(requireEmmy);

var _COLON_pi = pi;
var __COLON_pi = _pi;

var specialchars =
[["_SYMBOL_$1", /\'(\w+)/g],
["$1_DOT_$2", /(\d+)\.(\d+)/g],
["__COLON_$1", /\-:(\w+)/g],
["_COLON_$1", /:(\w+)/g],
["$1__GT_$2", /(\w+)->(\w+)/g],
["$1_GT__$2", /(\w+)<-(\w+)/g], //part3 added
[" state__GT_t ", /(\s+)time(\s+)/g],
["[ / -$1 $2 ],", /-(\w+)\/(\d+)/g],
["[ / $1 $2 ],", /(\w+)\/(\d+)/g]];

var mathfns =
[["_SLASH_", /\//g],
["_STAR_", /\*/g],
["_", /\-/g],
["_PLUS_", /\+/g],
["expt", /\^/g]];

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

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

requireEmmyRename("D","drv");

var specialchars_2 =
    [["_SYMBOL_$1^$2", /_SYMBOL_(\w+)expt(\w+)/],
     ["$1.$2", /(\d+)_DOT_(\d+)/g],
     ["-$1", /_(\d+)/g]];

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

["__GT_infix", "simplify"].map(requireEmmy);

var preAmble = (j) =>
j[0] != "define"
? ["__GT_infix", ["simplify", 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(edgeCases(letFlat(replaceD(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, 8) == "_SYMBOL_"
? 'symbol("' + j.substring(8, j.length) + '")'
:j.constructor == Array
? jsonToJs(j[0]) + ".call(null,"  +  j.slice(1).map(jsonToJs) + ")"
: j;

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

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

var show_expression = identity;
var velocities = velocity;
var coordinates = coordinate;
var vector_length = count; null;
    

Warm up

(+ 1 2 3) 
    

(cos :pi) 
    

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

3 Hamiltonian Mechanics

3.1 Hamilton’s Equations

Computing Hamilton's equations

(define ((Hamilton-equations Hamiltonian) q p)
(let ((state-path (qp->H-state-path q p)))
(- (D state-path)
(compose (Hamiltonian->state-derivative Hamiltonian)
state-path)))) 
    

(define ((Hamiltonian->state-derivative Hamiltonian) H-state)
(up 1
(((partial 2) Hamiltonian) H-state)
(- (((partial 1) Hamiltonian) H-state)))) 
    

(define ((qp->H-state-path q p) t)
(up t (q t) (p t))) 
    

(define ((H-rectangular m V) state)
(let ((q (coordinate state))
(p (momentum state)))
(+ (/ (square p) (* 2 m))
(V (ref q 0) (ref q 1))))) 
    

(show-expression
(let ((V (literal-function 'V (-> (X Real Real) Real)))
(q (up (literal-function 'x)
(literal-function 'y)))
(p (down (literal-function 'p_x)
(literal-function 'p_y))))
(((Hamilton-equations (H-rectangular 'm V)) q p) 't)))
    

(define (Legendre-transform F)
(let ((w-of-v (D F)))
(define (G w)
(let ((zero (compatible-zero w)))
(let ((M ((D w-of-v) zero))
(b (w-of-v zero)))
(let ((v (solve-linear-left M (- w b))))
(- (* w v) (F v))))))
G))
    

(define ((Lagrangian->Hamiltonian Lagrangian) H-state)
(let ((t (time H-state))
(q (coordinate H-state))
(p (momentum H-state)))
(define (L qdot)
(Lagrangian (up t q qdot)))
((Legendre-transform L) p)))
    

(define ((Hamiltonian->Lagrangian Hamiltonian) L-state)
(let ((t (time L-state))
(q (coordinate L-state))
(qdot (velocity L-state)))
(define (H p)
(Hamiltonian (up t q p)))
((Legendre-transform H) qdot)))
    

(define ((L-rectangular m V) local)
(let ((q (coordinate local))
(qdot (velocity local)))
(- (* 1/2 m (square qdot))
(V (ref q 0) (ref q 1)))))
    

(show-expression
((Lagrangian->Hamiltonian
(L-rectangular
’m
(literal-function ’V (-> (X Real Real) Real))))
(up ’t (up ’x ’y) (down ’p_x ’p_y))))
    

(define F
(literal-function ’F
(-> (UP Real (UP Real Real) (DOWN Real Real)) Real)))
    

(define G
(literal-function ’G
(-> (UP Real (UP Real Real) (DOWN Real Real)) Real)))
    

(define H
(literal-function ’H
(-> (UP Real (UP Real Real) (DOWN Real Real)) Real)))
    

((+ (Poisson-bracket F (Poisson-bracket G H))
(Poisson-bracket G (Poisson-bracket H F))
(Poisson-bracket H (Poisson-bracket F G)))
(up ’t (up ’x ’y) (down ’px ’py)))
    

(show-expression
((Lagrangian->Hamiltonian
(L-central-polar ’m (literal-function ’V)))
(up ’t (up ’r ’phi) (down ’p_r ’p_phi))))
    

(show-expression
(((Hamilton-equations
(Lagrangian->Hamiltonian
(L-central-polar ’m (literal-function ’V))))
(up (literal-function ’r)
(literal-function ’phi))
(down (literal-function ’p_r)
(literal-function ’p_phi)))
’t))
    

(define ((L-axisymmetric-top A C gMR) local)
(let ((q (coordinate local))
(qdot (velocity local)))
(let ((theta (ref q 0))
(thetadot (ref qdot 0))
(phidot (ref qdot 1))
(psidot (ref qdot 2)))
(+ (* 1/2 A
(+ (square thetadot)
(square (* phidot (sin theta)))))
(* 1/2 C
(square (+ psidot (* phidot (cos theta)))))
(* -1 gMR (cos theta))))))
    

(show-expression
((Lagrangian->Hamiltonian (L-axisymmetric-top ’A ’C ’gMR))
(up ’t
(up ’theta ’phi ’psi)
(down ’p_theta ’p_phi ’p_psi))))
    

(show-expression
((Lagrangian->Hamiltonian
(L-periodically-driven-pendulum ’m ’l ’g ’a ’omega))
(up ’t ’theta ’p_theta)))
    

(define (H-pend-sysder m l g a omega)
(Hamiltonian->state-derivative
(Lagrangian->Hamiltonian
(L-periodically-driven-pendulum m l g a omega))))
    

(define window (frame :-pi :pi -10.0 10.0))
    

(define ((monitor-p-theta win) state)
(let ((q ((principal-value :pi) (coordinate state)))
(p (momentum state)))
(plot-point win q p)))
    

(let ((m 1.0)	;m=1kg
(l 1.0)

;l=1m
(g 9.8)

;g=9.8m/s2
(A 0.1)

;A=1/10 m
(omega (* 2 (sqrt 9.8))))

((evolve H-pend-sysder m l g A omega)

(up 0.0

;t0=0
1.0

;theta0=1 rad
0.0)

;p0=0 kg m2/s
(monitor-p-theta window)

0.01

;plot interval
100.0

;final time
1.0e-12))
    

(define (driven-pendulum-map m l g A omega)
(let ((advance (state-advancer H-pend-sysder m l g A omega))
(map-period (/ :2pi omega)))
(lambda (theta ptheta return fail)
(let ((ns (advance
(up 0 theta ptheta)
; initial state
map-period)))
; integration interval
(return ((principal-value :pi)
(coordinate ns))
(momentum ns))))))
    

(define win (frame :-pi :pi -20 20))
    

(let ((m 1.0)
;m=1kg
(l 1.0)
;l=1m
(g 9.8)
;g=9.8m/s2
(A 0.05))
;A=1/20m
(let ((omega0 (sqrt (/ g l))))
(let ((omega (* 4.2 omega0)))
(explore-map
win
(driven-pendulum-map m l g A omega)
1000)))) ;1000 points for each initial condition
    

(define (HHmap E dt sec-eps int-eps)
(define ((make-advance advancer eps) s dt)
(advancer s dt eps))
(let ((adv
(make-advance (state-advancer HHsysder) int-eps)))
(lambda (y py cont fail)
(let ((initial-state (section->state E y py)))
(if (not initial-state)
(fail)
(find-next-crossing initial-state adv dt sec-eps
(lambda (crossing-state running-state)
(cont (ref (coordinate crossing-state)
1)
(ref (momentum crossing-state)
1)))))))))
    

(define (section->state E y py)
(let ((d (- E (+ (HHpotential (up 0 (up 0 y)))
(* 1/2 (square py))))))
(if (>= d 0.0)
(let ((px (sqrt (* 2 d))))
(up 0 (up 0 y) (down px py)))
#f)))
    

(define (HHHam s)
(+ (* 1/2 (square (momentum s)))
(HHpotential s)))
    

(define (HHpotential s)
(let ((x (ref (coordinate s) 0))
(y (ref (coordinate s) 1)))
(+ (* 1/2 (+ (square x) (square y)))
(- (* (square x) y) (* 1/3 (cube y))))))
    

(define (HHsysder)
(Hamiltonian->state-derivative HHHam))
    

(define (find-next-crossing state advance dt sec-eps cont)
(let lp ((s state))
(let ((next-state (advance s dt)))
(if (and (> (ref (coordinate next-state) 0) 0)
(< (ref (coordinate s) 0) 0))
(let ((crossing-state
(refine-crossing sec-eps advance s)))
(cont crossing-state next-state))
(lp next-state)))))
    

(define (refine-crossing sec-eps advance state)
(let lp ((state state))
(let ((x (ref (coordinate state) 0))
(xd (ref (momentum state) 0)))
(let ((zstate (advance state (- (/ x xd)))))
(if (< (abs (ref (coordinate zstate) 0))
sec-eps)
zstate
(lp zstate))))))
    

(define ((standard-map K) theta I return failure)
(let ((nI (+ I (* K (sin theta)))))
(return ((principal-value :2pi) (+ theta nI))
((principal-value :2pi) nI))))
    

(define window (frame 0.0 :2pi 0.0 :2pi))
(explore-map window (standard-map 0.6) 2000)