Solutions in search of problems

Suppose you have a function like (define foo (lambda (x) (- (* x x x) 30))) and you want to find x such that (foo x) = 0. There are a few ways to go about this. If you can find two different x such that (foo x) is positive for one and negative for the other, then (foo x) must be zero somewhere in between. A simple binary search will find it.

(define (bisection-method f left right)
  (let* ((midpoint (average left right))
         (fmid     (f midpoint)))
    (if (< (abs fmid) 1e-8)
        midpoint
        (let ((fl (f left))
              (fr (f right)))
          (cond ((same-sign? fl fr) (error "Left and right not on opposite sides."))
                ((same-sign? fmid fr) (bisection-method f left midpoint))
                ((same-sign? fl fmid) (bisection-method f midpoint right))
                (else (error "shouldn't happen")))))))

(define (average l r) (/ (+ l r) 2))

(define (same-sign? l r)
  (or (and (positive? l)
           (positive? r))
      (and (negative? l)
           (negative? r))))

1 ]=> (cos 2)

;Value: -.4161468365471424

1 ]=> (cos 1)

;Value: .5403023058681398

1 ]=> (bisection-method cos 1.0 2.0)
1. 2.
1.5 2.
1.5 1.75
1.5 1.625
1.5625 1.625
1.5625 1.59375
1.5625 1.578125
1.5703125 1.578125
1.5703125 1.57421875
1.5703125 1.572265625
1.5703125 1.5712890625
1.5703125 1.57080078125
1.570556640625 1.57080078125
1.5706787109375 1.57080078125
1.57073974609375 1.57080078125
1.570770263671875 1.57080078125
1.5707855224609375 1.57080078125
1.5707931518554687 1.57080078125
1.5707931518554687 1.5707969665527344
1.5707950592041016 1.5707969665527344
1.570796012878418 1.5707969665527344
1.570796012878418 1.5707964897155762
1.570796251296997 1.5707964897155762
1.570796251296997 1.5707963705062866
1.5707963109016418 1.5707963705062866
1.5707963109016418 1.5707963407039642
;Value: 1.570796325802803

Rather than selecting the midpoint between the two prior guesses, you can pretend that your function is linear between the guesses and interpolate where the zero should be. This can converge quicker.

(define (secant-method f x1 x2)
  (display x1) (display " ") (display x2) (newline)
  (let ((f1 (f x1))
        (f2 (f x2)))
    (if (< (abs f1) 1e-8)
        x1
        (let ((x0 (/ (- (* x2 f1) (* x1 f2))
                     (- f1 f2))))
          (secant-method f x0 x1)))))

1 ]=> (secant-method cos 0.0 4.0)
0. 4.
2.418900874126076 0.
1.38220688493168 2.418900874126076
1.5895160570280047 1.38220688493168
1.5706960159120333 1.5895160570280047
1.5707963326223677 1.5706960159120333
;Value: 1.5707963326223677

If you know the derivative of f, then you can use Newton's method to find the solution.

(define (newtons-method f f-prime guess)
  (display guess) (display " ") (newline)
  (let ((dy (f guess)))
    (if (< (abs dy) 1e-8)
        guess
        (let ((dy/dx (f-prime guess)))
          (newtons-method f f-prime (- guess (/ dy dy/dx)))))))

1 ]=> (newtons-method cos (lambda (x) (- (sin x))) 2.0)
2. 
1.5423424456397141 
1.5708040082580965 
1.5707963267948966 
;Value: 1.5707963267948966

Here's another example. We'll find the cube root of 30 by solving (lambda (x) (- (* x x x) 30)).

(define (cube-minus-thirty x) (- (* x x x) 30))

1 ]=> (bisection-method cube-minus-thirty 0.0 4.0)
0. 4.
2. 4.
3. 4.
3. 3.5
3. 3.25
3. 3.125
3.0625 3.125
3.09375 3.125
3.09375 3.109375
3.1015625 3.109375
3.10546875 3.109375
3.10546875 3.107421875
3.1064453125 3.107421875
3.10693359375 3.107421875
3.107177734375 3.107421875
3.107177734375 3.1072998046875
3.107177734375 3.10723876953125
3.107208251953125 3.10723876953125
3.1072235107421875 3.10723876953125
3.1072311401367187 3.10723876953125
3.1072311401367187 3.1072349548339844
3.1072311401367187 3.1072330474853516
3.107232093811035 3.1072330474853516
3.107232093811035 3.1072325706481934
3.1072323322296143 3.1072325706481934
3.107232451438904 3.1072325706481934
3.107232451438904 3.1072325110435486
3.107232481241226 3.1072325110435486
3.1072324961423874 3.1072325110435486
3.107232503592968 3.1072325110435486
3.107232503592968 3.1072325073182583
3.107232505455613 3.1072325073182583
3.107232505455613 3.1072325063869357
;Value: 3.1072325059212744

1 ]=> (secant-method cube-minus-thirty 0.0 4.0)
0. 4.
1.875 0.
8.533333333333333 1.875
2.1285182547245376 8.533333333333333
2.341649751209406 2.1285182547245376
3.4857887202177547 2.341649751209406
3.0068542655016235 3.4857887202177547
3.0957153766467633 3.0068542655016235
3.1076136741672546 3.0957153766467633
3.1072310897513415 3.1076136741672546
3.1072325057801455 3.1072310897513415
;Value: 3.1072325057801455

1 ]=> (define (cube-minus-thirty-prime x) (* 3 x x))

1 ]=> (newtons-method cube-minus-thirty cube-minus-thirty-prime 4.0)
4. 
3.2916666666666665 
3.1173734622300557 
3.10726545916981 
3.1072325063033337 
3.107232505953859 
;Value: 3.107232505953859