Skip to content

Instantly share code, notes, and snippets.

@ddrone
Created October 23, 2012 14:42
Show Gist options
  • Save ddrone/3939166 to your computer and use it in GitHub Desktop.
Save ddrone/3939166 to your computer and use it in GitHub Desktop.
#lang racket
(require plot)
(plot-new-window? #t)
(define (flatmap f ls)
(if (null? ls) '()
(append (f (car ls)) (flatmap f (cdr ls)))))
(define (range n)
(define (range-iter k)
(if (= k n) '()
(cons k (range-iter (+ 1 k)))))
(range-iter 0))
(define (float-range start step n)
(map (lambda (k)
(+ start (* step k)))
(range n)))
(define ((f r) x)
(* r x (- 1 x)))
(define (iter f n z)
(if (= n 0) z
(iter f (- n 1) (f z))))
(define (scanl f n z)
(if (= n 0) (list z)
(cons z (scanl f (- n 1) (f z)))))
(define rs (float-range 1 0.01 1000))
(define (get-points r)
(map (lambda (x) (vector r x)) (scanl (f r) 100 (iter (f r) 400 (random)))))
(define (show-plot)
(plot (points (flatmap get-points rs) #:sym 'dot)))
#lang racket
(require plot)
(define (flatmap f ls)
(if (null? ls) '()
(append (f (car ls)) (flatmap f (cdr ls)))))
(define (map-prod f ls1 ls2)
(flatmap (lambda (x) (map (lambda (y) (f x y)) ls2)) ls1))
(define (range n)
(define (range-iter k)
(if (= k n) '()
(cons k (range-iter (+ 1 k)))))
(range-iter 0))
(define (float-range start step n)
(map (lambda (k)
(+ start (* step k)))
(range n)))
(define axes-range (float-range -1.5 0.05 60)) ; 0.05 60 ; 0.02 150
(define test-points (map-prod make-rectangular axes-range axes-range))
(define roots (list (make-polar 1 0)
(make-polar 1 (* 2/3 pi))
(make-polar 1 (* 4/3 pi))))
(define (enumerate ls)
(define (enumerate-helper ls n)
(if (null? ls) '()
(cons (cons n (car ls))
(enumerate-helper (cdr ls) (+ 1 n)))))
(enumerate-helper ls 0))
(define (root-num z)
(min-pair (enumerate (map (lambda (r) (magnitude (- r z))) roots))))
(define (min-pair ls)
(cond
((null? ls) (error "list should not be empty"))
((null? (cdr ls)) (car ls))
(else (let ([rest-min (min-pair (cdr ls))])
(if (< (cdar ls) (cdr rest-min))
(car ls)
rest-min)))))
(define (f z)
(- (* z z z) 1))
(define (df z)
(* 3 z z))
(define (newton-iter z)
(- z (/ (f z) (df z))))
(define (iter f n z)
(if (= n 0) z
(iter f (- n 1) (f z))))
(define (scanl f n z)
(if (= n 0) (list z)
(cons z (scanl f (- n 1) (f z)))))
(define results
(map (lambda (x) (root-num (iter newton-iter 100 x))) test-points))
(define (complex->vector c)
(vector (real-part c)
(imag-part c)))
(define (extract n xs ys)
(cond
((null? ys) '())
((= n (caar ys)) (cons (complex->vector (car xs)) (extract n (cdr xs) (cdr ys))))
(else (extract n (cdr xs) (cdr ys)))))
(define class0 (extract 0 test-points results))
(define class1 (extract 1 test-points results))
(define class2 (extract 2 test-points results))
(define (particle-path z)
(map complex->vector (interpolate (scanl newton-iter 10 z))))
(define (interpolate ls)
(cond
((null? ls) ls)
((null? (cdr ls)) ls)
(else (append (float-range (car ls) (* 0.01 (- (cadr ls) (car ls))) 100)
(interpolate (cdr ls))))))
(define (show-plot)
(plot (list (points class0 #:sym 'fullcircle1 #:color 'red)
(points class1 #:sym 'fullcircle1 #:color 'green)
(points class2 #:sym 'fullcircle1 #:color 'blue)
(points (particle-path 1.0+1.2i) #:sym 'fullcircle1)
(points (particle-path -1.0-1.4i) #:sym 'fullcircle1))))
(show-plot)
#lang racket
(require plot)
(define (flatmap f ls)
(if (null? ls) '()
(append (f (car ls)) (flatmap f (cdr ls)))))
(define (range n)
(define (range-iter k)
(if (= k n) '()
(cons k (range-iter (+ 1 k)))))
(range-iter 0))
(define (float-range start step n)
(map (lambda (k)
(+ start (* step k)))
(range n)))
(define (float-range2 start stop n)
(float-range start (/ (- stop start) n) n))
(define ((f r) x)
(* r x (- 1 x)))
(define (id x) x)
(define (scanl f n z)
(if (= n 0) (list z)
(cons z (scanl f (- n 1) (f z)))))
(define ((show-path r) x)
(append (map (lambda (y) (vector x y)) (float-range2 x ((f r) x) 100))
(map (lambda (x1) (vector x1 ((f r) x))) (float-range2 x ((f r) x) 100))))
(define (particle-path r x)
(flatmap (show-path r) (scanl (f r) 100 x)))
(define (interpolate ls)
(cond
((null? ls) ls)
((null? (cdr ls)) ls)
(else (append (float-range (car ls) (* 0.01 (- (cadr ls) (car ls))) 100)
(interpolate (cdr ls))))))
(define (show-plot x0 r lx)
(plot (list (function (f r) 0 lx)
(function id 0 lx)
(points (particle-path r x0) #:sym 'fullcircle1))))
(show-plot 0.5 3.9 2)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment