SICP Exercise 2.14: Broken Intervals

After considerable work, Alyssa P. Hacker delivers her finished system. Several years later, after she has forgotten all about it, she gets a frenzied call from an irate user, Lem E. Tweakit. It seems that Lem has noticed that the formula for parallel resistors can be written in two algebraically equivalent ways:

  R1R2 
R1 + R2

and

     1     
1/R1 + 1/R2

He has written the following two programs, each of which computes the parallel-resistors formula differently:

(define (par1 r1 r2)
  (div-interval (mul-interval r1 r2)
                (add-interval r1 r2)))
(define (par2 r1 r2)
  (let ((one (make-interval 1 1))) 
    (div-interval one
                  (add-interval (div-interval one r1)
                                (div-interval one r2)))))

Lem complains that Alyssa's program gives different answers for the two ways of computing. This is a serious complaint.

Demonstrate that Lem is right. Investigate the behavior of the system on a variety of arithmetic expressions. Make some intervals A and B, and use them in computing the expressions A/A and A/B. You will get the most insight by using intervals whose width is a small percentage of the center value. Examine the results of the computation in center-percent form (see exercise 2.12).

Before we start building intervals and trying things out, let's first define a procedure for displaying intervals in center-percent form. The implementation I defined uses the center and percent selectors from exercise 2.12:

(define (print-center-percent i)
  (newline)
  (display (center i))
  (display " ± ")
  (display (percent i))
  (display "%"))

Now we can define a couple of intervals, a and b, and check this works okay:

> (define a (make-center-percent 100 1))
> (define b (make-center-percent 150 2))
> a
'(99 . 101)
> b
'(147 . 153)
> (print-center-percent a)
100 ± 1%
> (print-center-percent b)
150 ± 2%

All seems to be in order here, so let's try out par1 and par2:

> (print-center-percent (par1 a a))
50.02000200020002 ± 2.999200239928031%
> (print-center-percent (par2 a a))
50.0 ± 1.000000000000007%
> (print-center-percent (par1 a b))
60.05617438064145 ± 4.597193908142474%
> (print-center-percent (par2 a b))
59.99855963126561 ± 1.400072020165649%
> (print-center-percent (par1 b b))
75.12004801920769 ± 5.993607670795042%
> (print-center-percent (par2 b b))
75.0 ± 2.0%

So Lem is right - two equivalent formulae produce different results!

Note that if we change the tolerances to 0% for both values (i.e. we just use the center values for calculating the resistance for the parallel circuits), we find that R_aa = 50Ω, R_ab = 60Ω and R_bb = 75Ω. So not only do they produce different results, but par2 both produces results with center values closer to the values we might expect and produces tighter tolerances.

Now onto A/A and A/B. Let's also throw B/A and B/B into the mix as well:

> (print-center-percent (div-interval a a))
1.0002000200020003 ± 1.9998000199980077%
> (print-center-percent (div-interval a b))
0.6670668267306923 ± 2.999400119976001%
> (print-center-percent (div-interval b a))
1.5004500450045004 ± 2.9994001199760123%
> (print-center-percent (div-interval b b))
1.0008003201280513 ± 3.9984006397441028%

The first point of interest to note is that our current system has no idea of identity. If a represents a particular instance of an interval then, while we may be unsure of its exact value, it does have an exact value and so (div-interval a a) should evaluate to 1 ± 0%. However it doesn't treat them as identical instances, but instead as two separate intervals that happen to have the same center value and tolerance, and so this doesn't happen.

The second point of interest is, with small tolerances at least, dividing one interval by another produces a new interval with a center point that is approximately equal to (actually slightly higher than) the value produced by dividing the center value of the numerator interval by the center value of the denominator interval and with a tolerance approximately equal to (actually slightly lower than) the sum of the tolerances of the two original intervals.

Assuming we use c_i to represent the center value of an interval and t_i to represent the percentage tolerance we can express this more succinctly as:

cx/y ≈ cx/cy ∧ cx/y ≥ cx/cy
tx/y ≈ tx/ty ∧ tx/y ≤ tx + ty

What about some other operations:

> (print-center-percent (mul-interval a a))
10001 ± 20000/10001%
> (print-center-percent (mul-interval a b))
15003 ± 5000/1667%
> (print-center-percent (mul-interval b a))
15003 ± 5000/1667%
> (print-center-percent (mul-interval b b))
22509 ± 10000/2501%
>
> (print-center-percent (add-interval a a))
200 ± 1%
> (print-center-percent (add-interval a b))
250 ± 8/5%
> (print-center-percent (add-interval b a))
250 ± 8/5%
> (print-center-percent (add-interval b b))
300 ± 2%
>
> (print-center-percent (sub-interval a a))
0 ± 1/0%
> (print-center-percent (sub-interval a b))
-50 ± -8%
> (print-center-percent (sub-interval b a))
50 ± 8%
> (print-center-percent (sub-interval b b))
0 ± 1/0%

I'll translate some of those percentages...

• 20000/10001% ≈ 1.9998%
• 5000/1667% ≈ 2.9994%
• 10000/2501% ≈ 3.9984%
• 8/5% = 1.6%

So we can see some further relationships here:

cxy ≈ cxcy ∧ cxy ≥ cxcy
txy ≈ txty ∧ txy ≤ tx + ty

cx+y = cx + cy
tx+y = (tx + ty) / 2

cx-y = cx - cy
tx-y = 100(wx + wy)/(cx - cy)

Okay, so I "cheated" with t_x-y and actually worked out the formula from the definitions of centers, widths tolerance and sub-interval. However hopefully all of the above gives us some indication as to what's going on with Lem's programs. While the two programs are algebraically equivalent, the calculations they perform are different. As a result the center values and tolerances of the intervals they produce will also be different. In fact this is perhaps hinting at the answer to the next exercise...