source: project/release/4/derivative/trunk/derivative.scm @ 30617

Last change on this file since 30617 was 30617, checked in by Ivan Raikov, 7 years ago

initial import of derivative, a library for computing the symbolic derivative of an arithmetic expression

File size: 5.7 KB
Line 
1
2(module derivative
3       
4        (expr? fold-expr eval-expr derivative-expr
5         newton-search )
6
7        (import scheme chicken)
8        (require-extension matchable)
9
10;; Expression type
11
12(define (binop? x)
13  (case x 
14    ((+ - * / pow) #t)
15    (else #f)))
16
17
18(define (unop? x)
19  (case x 
20    ((sin cos ln neg exp deriv const) #t)
21    (else #f)))
22
23
24(define (expr? x)
25  (match x 
26         ((? number?) #t)
27         ((? symbol?) #t)
28         (((? binop?) (? expr?) (? expr?)) #t)
29         (((? unop?) (? expr?)) #t)
30         (else #f)
31         ))
32
33
34(define (fold-expr fatom funary fbinary expr)
35  (match expr
36         ((or (? symbol?) (? number?)) 
37          (fatom expr))
38         (((and op (? binop?)) (and (? expr?) l) (and (? expr?) r))
39          (fbinary op (fold-expr fatom funary fbinary l)
40                   (fold-expr fatom funary fbinary r)))
41         (((and op (? unop?)) (and (? expr?) e))
42          (funary op (fold-expr fatom funary fbinary e)))
43         ))
44
45
46;; Does the given expression contains a variable (symbol)
47
48(define (has-var? expr)
49  (fold-expr symbol? (lambda (op x) x) (lambda (op x y) (or x y)) expr))
50
51;; Evaluates the given expression with the given input.
52
53(define (eval-expr var e)
54  (define (eval-binop op left right)
55    (case op 
56      ((+) (+ left right))
57      ((-) (- left right))
58      ((*) (* left right))
59      ((/) (/ left right))
60      ((pow) (expt left right))
61      ))
62  (define (eval-unop op e)
63    (case op 
64      ((sin) (sin e))
65      ((cos) (cos e))
66      ((ln)  (log e))
67      ((neg) (- e))
68      ((exp) (exp e))
69      ((deriv) 1.0)
70      ((const) (error 'eval-expr "cannot evaluate symbolic constant" e))
71      ))
72  (define (eval-atom x)
73    (cond ((number? x) x)
74          ((symbol? x) var)
75          (else (error 'eval-expr "invalid atom" x))
76          ))
77  (fold-expr eval-atom eval-unop eval-binop e)
78  )
79
80;; TODO:
81;;    ~~f = f
82;;    ~f + g = g - f
83;;    f + (~g) = f - g
84;;    f - (~g) = f + g
85;;    ~f - g = ~(f + g)
86;;    (~f) * g = ~(f * g)
87;;    f * (~g) = ~(f * g)
88;;    (f/g) * h = (f * h) / g
89;;    f * (g/h) = (f * g) / h
90;;    (~f) / g = ~(f / g)
91;;    f / (~g) = ~(f / g)
92;;    (f / g) / h = f / (g * h)
93;;    f / (g / h) = (f * h) / g
94;;    f ^ (~g) = 1 / (f ^ g)
95;;    (f ^ g) ^ h = f ^ (g * h)
96
97
98(define (simplify-expr e)
99  (define (simplify-binop op left right)
100    (match (list op left right)
101      (('+ (and (? number?) x) (and (? number?) y)) (+ x y))
102      (('+ 0.0 y) y)
103      (('+ x 0.0) x)
104      (('- (and (? number?) x) (and (? number?) y)) (- x y))
105      (('- 0.0 y) `(neg ,y))
106      (('- x 0.0) x)
107      (('* (and (? number?) x) (and (? number?) y)) (* x y))
108      (('* (and (? number?) x) ('* (and (? number?) y) rst)) 
109       `(* ,(* x y) ,rst))
110      (('* x 1.0) x)
111      (('* 1.0 x) x)
112      (('* 0.0 x) 0.0)
113      (('* x 0.0) 0.0)
114      (('/ (and (? number?) x) (and (? number?) y)) (/ x y))
115      (('/ x 1.0) x)
116      (('pow x 1.0) x)
117      (('cos 0.0) 1.0)
118      (('sin 0.0) 0.0)
119      (('ln 1.0)  0.0)
120      (else `(,op ,left ,right))
121      ))
122  (define (simplify-unop op e)
123    (case op 
124      ((neg) (if (number? e) (- e) `(,op ,e)))
125      (else `(,op ,e))
126      ))
127  (define (simplify-atom x)
128    (cond ((number? x) x)
129          ((symbol? x) x)
130          (else (error 'simplify-expr "invalid atom" x))
131          ))
132  (fold-expr simplify-atom simplify-unop simplify-binop e)
133  )
134
135
136;; Computes the derivative of the given expression.
137
138(define (derivative-expr e)
139  (define (fatom x)
140    (cond ((number? x) (cons x 0.0))
141          ((symbol? x) (cons x `(deriv ,x)))
142          (else (error 'derivative-expr "invalid atom" e)))
143    )
144  (define (fbinary op ldl rdr)
145    (match-let (((l . dl) ldl)
146                ((r . dr) rdr))
147       (case op
148         ((pow)
149          (let ((expr  `(pow ,l ,r)))
150            (if (has-var? r)
151                (let* ((frac  `(/ (* ,r ,dl) ,l))
152                       (inner `(+ (* ,dr (ln ,l)) ,frac)))
153                  (cons expr `(* ,expr ,inner)))
154                (let ((new `(- ,r 1.0)))
155                  (cons expr `(* ,r (* (pow ,l ,new) ,dl))))
156                ))
157          )
158         ((+)  (cons `(+ ,l ,r) `(+ ,dl ,dr)))
159         ((-)  (cons `(- ,l ,r) `(- ,dl ,dr)))
160         ((*)  (cons `(* ,l ,r) `(+ (* ,l ,dr) (* ,dl ,r))))
161         ((/)  (cons `(/ ,l ,r) 
162                     `(/ (- (* ,dl ,r) (* ,l ,dr)) (pow ,r 2.0))))
163         )))
164  (define (funary op xdx)
165    (match-let (((x . dx) xdx))
166       (case op 
167         ((sin) (cons `(sin ,x) `(* (cos ,x) ,dx)))
168         ((cos) (cons `(cos ,x) `(* (neg (sin ,x)) ,dx)))
169         ((ln)  (cons `(ln ,x)  `(/ ,dx ,x)))
170         ((neg) (cons `(neg ,x) `(neg ,dx)))
171         ((exp) (cons `(exp ,x) `(* ,dx (exp ,x))))
172         ((const) (cons x 0.0))
173         )))
174  (simplify-expr (cdr (fold-expr fatom funary fbinary e)))
175  )
176
177
178(define machine-epsilon
179  (let loop ((e 1.0))
180     (if (= 1.0 (+ e 1.0))
181         (* 2 e)
182         (loop (/ e 2)))))
183
184
185(define (close-enough? h1 h2 tolerance)
186  (<= (magnitude (- h1 h2))
187      (* .5 (max tolerance machine-epsilon)
188         (+ (magnitude h1) (magnitude h2) 2.0))))
189
190
191;; Newton’s method for finding the root of a function f(y): if
192;; y0 is an approximation to the root of f, then:
193;;
194;; y0 - f(y0)/f'(y0)
195;;
196;; is a better approximation, where f'(y0) is the first derivative of
197;; f evaluation at y0.
198;;
199;; Argument eps is the error tolerance, maxit is the maximum number of
200;; iterations allowed.
201;;
202(define (newton-search fdf x0 eps maxit) 
203  (define (newton-improve xn)
204    (let ((fn (eval-expr xn (car fdf)))
205          (dfn (eval-expr xn (cdr fdf))))
206      (- xn (/ fn dfn))))
207  (let lp ((xn x0) (it 0))
208    (let ((xn+1 (newton-improve xn)))
209      (if (close-enough? xn xn+1 eps)
210          (/ (+ xn xn+1) 2)
211          (if (< it maxit)
212              (lp xn+1 (+ it 1))
213              #f))
214      ))
215  )
216
217
218
219)
Note: See TracBrowser for help on using the repository browser.