aboutsummaryrefslogtreecommitdiffstats
path: root/work
diff options
context:
space:
mode:
Diffstat (limited to 'work')
-rw-r--r--work/funcdiff_play.scm212
1 files changed, 212 insertions, 0 deletions
diff --git a/work/funcdiff_play.scm b/work/funcdiff_play.scm
new file mode 100644
index 0000000..429933c
--- /dev/null
+++ b/work/funcdiff_play.scm
@@ -0,0 +1,212 @@
+
+;;; Jan 29th
+
+(load "~/thesis/scmutils/src/calculus/load.scm")
+
+;(define R2 (make-manifold R^n-type 2)) ; doesn't work, so...
+(define R2 (rectangular 2))
+;(define U (patch 'origin R2)) ; also undefined so going to AIM-2005-003.pdf
+
+(define R2 (rectangular 2))
+(define P2 (polar/cylindrical 2))
+
+
+(define R2-chi-inverse (R2 '->point))
+(define R2-chi (R2 '->coords))
+(define P2-chi (P2 '->coords))
+(define P2-chi-inverse (P2 '->point))
+
+;;; Feb 1st
+
+(print-expression
+ ((compose R2-chi-inverse P2-chi)
+ (up 'x0 'y0)))
+;= (up (sqrt (+ (expt x0 2) (expt y0 2))) (atan y0 x0))
+
+(print-expression
+ ((compose R2-chi P2-chi-inverse)
+ (up 'r0 'theta0)))
+
+; (up (* r0 (cos theta0)) (* r0 (sin theta0)))
+
+(define R2->R (-> (UP Real Real) Real))
+
+(define f
+ (compose (literal-function 'f-rect R2->R)
+ R2-chi))
+; there is a simpler syntax...
+
+(define R2-point (R2-chi-inverse (up 'x0 'y0)))
+(define P2-point (P2-chi-inverse (up 'r0 'theta0)))
+
+(print-expression (f R2-point))
+; (f-rect (up x0 y0))
+(print-expression (f P2-point))
+; (f-rect (up (* r0 (cos theta0)) (* r0 (sin theta0))))
+
+(define g (literal-manifold-function 'g-polar P2))
+
+(print-expression (g R2-point))
+; (g-polar (up (sqrt (+ (expt x0 2) (expt y0 2))) (atan y0 x0)))
+
+(instantiate-coordinates R2 '(x y))
+(instantiate-coordinates P2 '(r theta))
+
+(print-expression (x R2-point))
+;x0
+
+(print-expression (theta R2-point))
+;(atan y0 x0)
+
+(define h (+ (* x (square r)) (cube y)))
+
+(print-expression (h P2-point))
+; (+ (* (expt r0 3) (expt (sin theta0) 3)) (* (expt r0 3) (cos theta0)))
+
+(print-expression (h R2-point))
+; (+ (expt x0 3) (* x0 (expt y0 2)) (expt y0 3))
+
+(print-expression
+ (D h))
+; a-euclidean-derivative
+
+; pe is print-expression
+
+(pe ((D h) R2-point))
+; (down (+ (* 3 (expt x0 2)) (expt y0 2)) (+ (* 2 x0 y0) (* 3 (expt y0 2))))
+
+;(pe (D (h R2-point)))
+; ERROR
+
+(pe ((D (compose h R2-chi-inverse)) R2-point))
+; (down (+ (* 3 (expt x0 2)) (expt y0 2)) (+ (* 2 x0 y0) (* 3 (expt y0 2))))
+
+; TODO: formal definition of operator vs. function?
+
+;;; Feb 2
+
+(define (vector-field-procedure components coordinate-system)
+ (define (the-procedure f)
+ (compose (* (D (compose f (coordinate-system '->point)))
+ components)
+ (coordinate-system '->coords)))
+ the-procedure)
+
+(define (components->vecotr-field components coordinate-system)
+ (procedure->vector-field
+ (vector-field-procedure components coordinate-system)))
+
+(define v
+ (components->vector-field
+ (up (literal-function 'vx (-> (UP Real Real) Real))
+ (literal-function 'vy (-> (UP Real Real) Real)))
+ R2))
+
+; shorthand:
+; (define w (literal-vector-field 'v R2))
+
+;(pe ((v (literal-manifold-function 'f R2)) R2-point))
+
+(define (coordinatize vector-field coordsys)
+ (define ((v f) x)
+ (let ((b (compose (vector-field (coordsys '->coords))
+ (coordsys '->point))))
+ (* ((D f) x) (b x))))
+ (make-operator v))
+
+#|
+(pe (((coordinatize v R2)
+ (literal-function 'f (-> (UP Real Real) Real)))
+ (up 'x0 'y0)))
+|#
+
+; Note: nice summary of vector field properties on p12
+
+;(pe ((d/dx (square r)) R2-point))
+;(pe ((d/dx (square r)) P2-point))
+
+(define J (- (* x d/dy) (* y d/dx)))
+
+(series:for-each pe
+ (((exp (* 'a J)) R2-chi)
+ ((R2 '->point) (up 1 0)))
+ 6)
+;(up 1 0)
+;(up 0 a)
+;(up (* -1/2 (expt a 2)) 0)
+;(up 0 (* -1/6 (expt a 3)))
+;(up (* 1/24 (expt a 4)) 0)
+;(up 0 (* 1/120 (expt a 5)))
+
+; now do evolution on coordinates
+
+(define ((((evolution order)
+ delta-t vector-field)
+ manifold-function)
+ manifold-point)
+ (series:sum
+ (((exp (* delta-t vector-field))
+ manifold-function)
+ manifold-point)
+ order))
+
+(pe ((((evolution 6) 'a J) R2-chi)
+ ((R2 '->point) (up 1 0))))
+
+(pe ((((evolution 6) 2. J) R2-chi)
+ ((R2 '->point) (up 1 0))))
+
+#|
+
+(define mywindow (frame -4 4 -4 4))
+
+(define (plot-evolution win order initial a step)
+ (letrec
+ ((dostep
+ (lambda (val)
+ (cond
+ ((< val a) (let ((this-point ((((evolution order) val J) R2-chi)
+ ((R2 '->point) initial))))
+ (plot-point win
+ (time this-point)
+ (coordinate this-point))
+ (dostep (+ val step))))))))
+ (dostep 0.)))
+
+(plot-evolution mywindow 6 (up 1. 0.) 6. .01)
+
+
+
+(define (explore-evolution window order length)
+ (define (iterate-mything i x y)
+ (if (< i length)
+ (let ((this-point ((((evolution order) i J) R2-chi)
+ ((R2 '->point) (up x y)))))
+ (plot-point window (time this-point) (coordinate this-point))
+ (iterate-mything (+ .01 i) x y))
+ (button-loop x y)))
+ (define (button-loop ox oy)
+ (pointer-coordinates
+ window
+ (lambda (x y button)
+ (let ((temp button))
+ (cond ((eq? temp 0) (write-line (cons x (cons y (quote ()))))
+ (display " started.")
+ (iterate-mything 0 x y))
+ ((eq? temp 1) (write-line (cons ox (cons oy (quote ()))))
+ (display " continued.")
+ (iterate-mything 0 ox oy))
+ ((eq? temp 2) (write-line (cons x (cons y (quote ()))))
+ (display " hit.")
+ (button-loop ox oy)))))))
+ (newline)
+ (display "Left button starts a trajectory.")
+ (newline)
+ (display "Middle button continues a trajectory.")
+ (newline)
+ (display "Right button interrogates coordinates.")
+ (button-loop 0. 0.))
+
+(explore-evolution mywindow 5 .4)
+
+|#