#lang racket (require plot) ;;; ;;; This implements plotting a 2d graph of a one variable function. ;;; Singularities are handled properly (I think). ;;; The points evaluated are chosen using an adaptive strategy: ;;; - in areas where the graph is smooth, few points are chosen ;;; - in areas where the graph is oscillating, more points are chosen ;;; See http://goo.gl/Hdi9q (plot.lisp in the Maxima source) ;;; for code implementing the adaptive strategy. ;;; See chapeter 4.1 in the YACAS book of algorithms for ;;; an explanation of the algorithm: ;;; http://yacas.sourceforge.net/Algochapter4.html#c4s1 ; These count functions are used to measure the ; number of evaluations of the function to be drawn. (define count 0) (define (reset-count) (set! count 0)) (define (increase-count) (set! count (+ count 1))) (define (wrap f [value-returned-on-error #f]) (λ (x) (with-handlers ([(λ e #t) (λ x value-returned-on-error)]) (increase-count) (f x)))) (define (hill-or-valley? a b c) ; (x1,a) (x2,b) (x3,c) points with a b a) (> b c)))) (define (oscilates? a b c d e) ; (x1,a) (x2,b) (x3,c) (x4,d) (x5,e) ; points with a (list (list x ...) ...) ; The list xs is split into sublists. ; For all neighbors x1 and x2 in xs, (pred x1 x2) determines whether the list is split. ; Example: (split-between = (list 1 1 2 3 3)) => '((1 1) (2) (3 3)) (define (split-between pred xs) (let loop ([xs xs] [ys '()] [xss '()]) (match xs [(list) (reverse (cons (reverse ys) xss))] [(list x) (reverse (cons (reverse (cons x ys)) xss))] [(list x1 x2 more ...) (if (pred x1 x2) (loop more (list x2) (cons (reverse (cons x1 ys)) xss)) (loop (cons x2 more) (cons x1 ys) xss))]))) (define (plot2d unwrapped-f [x-min -5] [x-max 5] [y-min -5] [y-max 5]) ; wrap the function to be drawn, s.t. it ; returns #f in error situations (define f (wrap unwrapped-f #f)) (define (clipped-to-different-sides? p q) ; are the points p and q on different ; sides of the y-min / y-max limit ? (define py (vector-ref p 1)) (define qy (vector-ref q 1)) (or (not (number? py)) (not (number? qy)) (and (< py y-min) (> qy y-max)) (and (> py y-max) (< qy y-min)))) (define (remove-non-numbers ps) (filter (λ (p) (number? (vector-ref p 1))) ps)) ; 29 is a good value according to plot.lisp (define number-of-regions 29) ; region widh (define delta (/ (- x-max x-min) number-of-regions)) ; generate points by dividing the interval ; from x-min to x-max into number-of-regions regions, ; and calling region, which calls adaptive. (define points (append* (for/list ([i number-of-regions]) (define a (+ x-min (* delta i))) ; avoid rounding error and compute c as: (define c (+ x-min (* delta (+ i 1)))) (if (= i 0) ; keep the first point for the first region (region f a c) ; otherwise remove the first point (which is ; present as the end point if the preceding region) (rest (region f a c)))))) ; Split the point list into groups. Split between two points ; if they are on different sides of y-min and y-max. ; All connected points will be drawn with (lines ...). (define connected-points (split-between clipped-to-different-sides? points)) ; Display both the adaptive plot and the original plot ; for visual comparision. (displayln (list (begin0 (plot (list (map lines (map remove-non-numbers connected-points))) #:x-min x-min #:x-max x-max #:y-min y-min #:y-max y-max) (displayln (format "adaptive number of evaluations: ~a" count)) (reset-count)) (plot (function (wrap unwrapped-f +inf.0) x-min x-max) #:y-min y-min #:y-max y-max))) (displayln (format "original number of evaluations: ~a\n\n" count)) (reset-count)) ; Examples (plot2d /) ; check division by zero and singularity at x=0 (plot2d tan) ; check multiple singularities (plot2d (λ (x) (- (/ x)))) ; check symmetry (plot2d (λ (x) (/ (abs x)))) ; check connected components of same sign (plot2d (λ (x) (* 1e6 x))) ; check that lines with high slopes are drawn (plot2d (λ (x) (* 1e50 x))) ; check that lines with high slopes are drawn