Permalink
Browse files

curious features of interval arithmetic

  • Loading branch information...
1 parent 2aac943 commit cc66732eba73b472520e0a9521edb8da115d39a1 @sarabander committed Aug 4, 2011
Showing with 493 additions and 0 deletions.
  1. +56 −0 2.1/2.09.scm
  2. +20 −0 2.1/2.10.scm
  3. +107 −0 2.1/2.11.scm
  4. +23 −0 2.1/2.12.scm
  5. +54 −0 2.1/2.13.scm
  6. +52 −0 2.1/2.14.scm
  7. +47 −0 2.1/2.15.scm
  8. +134 −0 2.1/2.16.scm
View
56 2.1/2.09.scm
@@ -0,0 +1,56 @@
+
+(define (add-interval x y)
+ (make-interval (+ (lower-bound x) (lower-bound y))
+ (+ (upper-bound x) (upper-bound y))))
+
+(define (sub-interval x y)
+ (let ((p1 (- (lower-bound x) (lower-bound y)))
+ (p2 (- (lower-bound x) (upper-bound y)))
+ (p3 (- (upper-bound x) (lower-bound y)))
+ (p4 (- (upper-bound x) (upper-bound y))))
+ (make-interval (min p1 p2 p3 p4)
+ (max p1 p2 p3 p4))))
+
+(define (mul-interval x y)
+ (let ((p1 (* (lower-bound x) (lower-bound y)))
+ (p2 (* (lower-bound x) (upper-bound y)))
+ (p3 (* (upper-bound x) (lower-bound y)))
+ (p4 (* (upper-bound x) (upper-bound y))))
+ (make-interval (min p1 p2 p3 p4)
+ (max p1 p2 p3 p4))))
+
+(define (div-interval x y)
+ (mul-interval x
+ (make-interval (/ 1.0 (upper-bound y))
+ (/ 1.0 (lower-bound y)))))
+
+(define (width interval)
+ (/ (- (upper-bound interval)
+ (lower-bound interval))
+ 2.0))
+
+(define interval1 (make-interval 4.4 4.9))
+(define interval2 (make-interval 7.2 7.6))
+
+(define interval3 (make-interval 3.4 3.9))
+(define interval4 (make-interval 9.2 9.6))
+
+(width interval1)
+(width interval3) ; => both 0.25
+(width interval2)
+(width interval4) ; => both 0.2
+
+;; Width in addition and subtraction depends only on argument interval widths
+(width (add-interval interval1 interval2)) ; => 0.45
+(width (sub-interval interval1 interval2)) ; => 0.45
+
+(width (add-interval interval3 interval4)) ; => 0.45
+(width (sub-interval interval3 interval4)) ; => 0.45
+
+;; Width in multiplication and division depends also on interval boundaries
+(width (mul-interval interval1 interval2)) ; => 2.78
+(width (div-interval interval1 interval2)) ; => 0.051
+
+(width (mul-interval interval3 interval4)) ; => 3.08
+(width (div-interval interval3 interval4)) ; => 0.035
+
View
20 2.1/2.10.scm
@@ -0,0 +1,20 @@
+
+(define (div-interval x y)
+ (if (opposite-sign? (lower-bound y) (upper-bound y))
+ (printf "I won't divide by an interval spanning zero.\n")
+ (mul-interval x
+ (make-interval (/ 1.0 (upper-bound y))
+ (/ 1.0 (lower-bound y))))))
+
+(define (opposite-sign? a b)
+ (xor (negative? a) (negative? b)))
+
+(define (xor x y)
+ (or (and x (not y))
+ (and y (not x))))
+
+(div-interval (make-interval -8.3 -7.9)
+ (make-interval -2.3 1.4)) ; => error
+
+(div-interval (make-interval -8.3 -7.9)
+ (make-interval 1.1 1.4)) ; => '(-7.55 . -5.64)
View
107 2.1/2.11.scm
@@ -0,0 +1,107 @@
+
+;; Let x and y designate the two intervals we wish to multiply.
+;; We will use these prefixes:
+
+;; l - lower-bound (e.g. lx is the lower-bound of x)
+;; u - upper-bound
+
+;; Interval endpoints could be positive or negative. There are 16
+;; combinations when taking two intervals with four endpoints. Some of
+;; these are not "legal" in the sense that one or both of the intervals
+;; have their upper bound smaller than lower bound (+ − combination).
+;; We can immediately dismiss seven of such illegal cases. Remains nine.
+
+;; We pick the two multiplications to produce result's lower and upper bound
+;; by carefully considering, which of them gives the smallest and which the
+;; largest value. Only one combination of signs gives an ambiguous result,
+;; so we need four multiplications, and must use min and max functions.
+
+;; Next table illustrates the multiplication pattern.
+
+;;
+;; x y multiplicands
+;; ----- ----- ---------------
+;; legal? lx ux ly uy case lower upper
+;; ------ -- -- -- -- ---- ----- -----
+;;
+;; ╭───u────╮
+;; ✓ (−)(−) (−)(−) 1. [ux⋅uy, lx⋅ly]
+;; ╰────l───╯
+;;
+;; ╭───u────╮
+;; ✓ (−)(−) (−)(+) 2. [lx⋅uy, lx⋅ly]
+;; ╰─────l─────╯
+;;
+;; ✗ (−)(−) (+)(−)
+;;
+;; ╭──u──╮
+;; ✓ (−)(−) (+)(+) 3. [lx⋅uy, ux⋅ly]
+;; ╰─────l─────╯
+;;
+;; ╭───u────╮
+;; ✓ (−)(+) (−)(−) 4. [ux⋅ly, lx⋅ly]
+;; ╰──l──╯
+;;
+;; ╭────────┬──╮
+;; ✓ (−)(+) (−)(+) 5. 4 multiplications, choose min, max
+;; ╰─────┴──╯
+;;
+;; ✗ (−)(+) (+)(−)
+;;
+;; ╭───u────╮
+;; ✓ (−)(+) (+)(+) 6. [lx⋅uy, ux⋅uy]
+;; ╰─────l─────╯
+;;
+;; ✗ (+)(−) (−)(−)
+;; ✗ (+)(−) (−)(+)
+;; ✗ (+)(−) (+)(−)
+;; ✗ (+)(−) (+)(+)
+;;
+;; ╭─────u─────╮
+;; ✓ (+)(+) (−)(−) 7. [ux⋅ly, lx⋅uy]
+;; ╰──l──╯
+;;
+;; ╭────u───╮
+;; ✓ (+)(+) (−)(+) 8. [ux⋅ly, ux⋅uy]
+;; ╰──l──╯
+;;
+;; ✗ (+)(+) (+)(−)
+;;
+;; ╭────u───╮
+;; ✓ (+)(+) (+)(+) 9. [lx⋅ly, ux⋅uy]
+;; ╰───l────╯
+
+
+(define (make-interval a b)
+ (if (<= a b)
+ (cons a b)
+ (printf "Lower bound must be at most equal to upper bound.\n")))
+
+(define (mul-interval x y)
+ (let ((lx (lower-bound x)) (ux (upper-bound x))
+ (ly (lower-bound y)) (uy (upper-bound y))
+ (-? negative?) (+? positive?) (mint make-interval))
+ (cond ((and (-? lx) (-? ux)) ; cases:
+ (cond ((and (-? ly) (-? uy)) (mint (* ux uy) (* lx ly))) ; 1.
+ ((and (-? ly) (+? uy)) (mint (* lx uy) (* lx ly))) ; 2.
+ ((and (+? ly) (+? uy)) (mint (* lx uy) (* ux ly))))) ; 3.
+ ((and (-? lx) (+? ux))
+ (cond ((and (-? ly) (-? uy)) (mint (* ux ly) (* lx ly))) ; 4.
+ ((and (-? ly) (+? uy)) ; 5.
+ (let ((p1 (* lx ly)) (p2 (* lx uy))
+ (p3 (* ux ly)) (p4 (* ux uy)))
+ (mint (min p1 p2 p3 p4) (max p1 p2 p3 p4))))
+ ((and (+? ly) (+? uy)) (mint (* lx uy) (* ux uy))))) ; 6.
+ ((and (+? lx) (+? ux))
+ (cond ((and (-? ly) (-? uy)) (mint (* ux ly) (* lx uy))) ; 7.
+ ((and (-? ly) (+? uy)) (mint (* ux ly) (* ux uy))) ; 8.
+ ((and (+? ly) (+? uy)) (mint (* lx ly) (* ux uy)))))))) ; 9.
+
+(define i5 (make-interval 4.4 4.9))
+(define i6 (make-interval 7.2 7.6))
+
+(define i7 (make-interval 3.4 3.9))
+(define i8 (make-interval 9.2 9.6))
+
+(mul-interval i5 i6)
+(mul-interval i7 i8)
View
23 2.1/2.12.scm
@@ -0,0 +1,23 @@
+;; From the book
+(define (make-center-width c w)
+ (make-interval (- c w) (+ c w)))
+
+(define (center i)
+ (/ (+ (lower-bound i) (upper-bound i)) 2))
+
+(define (width i)
+ (/ (- (upper-bound i) (lower-bound i)) 2))
+;; ----------------------
+
+(define (make-center-percent c p)
+ (make-center-width c (* c p)))
+
+(define (percent i)
+ (/ (width i) (center i)))
+
+;; Tests
+(define 5R10% (make-center-percent 5.0 0.1))
+
+(percent 5R10%) ; => 0.1
+(width 5R10%) ; => 0.5
+(center 5R10%) ; => 5.0
View
54 2.1/2.13.scm
@@ -0,0 +1,54 @@
+
+;; Depends on the two preceding exercises
+
+;; Let x, y and z be intervals. We will use the following prefixes with them:
+
+;; l - lower-bound (e.g. lx is the lower-bound of x)
+;; u - upper-bound
+;; c - center of interval
+;; w - width of interval (distance from center to either bound)
+;; p - percentage tolerance
+;; w u - l u + l
+;; We need these formulas: p = ─ => w = p⋅c, w = ─────, c = ─────
+;; c 2 2
+;; u - l
+;; It follows that p = ─────, l = c - w and u = c + w
+;; u + l
+;;
+;; Let's construct a multiplication z = x⋅y. Assuming x > 0 and y > 0
+;; (see 2.11), and expanding x, y and z to interval notation, we see that
+;; z = x⋅y => [lz, uz] = [lx, ux]⋅[ly, uy] = [lx⋅ly, ux⋅uy], so
+;; lz = lx⋅ly and uz = ux⋅uy.
+;; uz - lz ux⋅uy - lx⋅ly
+;; What is pz in terms of px and py? Well, pz = ─────── = ─────────────.
+;; uz + lz ux⋅uy + lx⋅ly
+;;
+;; We shall express lx, ly, ux and uy in terms of cx, cy, px and py.
+;; Using the formulas, we get:
+;;
+;; lx = cx - wx = cx - px⋅cx = cx(1 - px)
+;; ly = cy - wy = cy - py⋅cy = cy(1 - py)
+;; ux = cx + wx = cx + px⋅cx = cx(1 + px)
+;; uy = cy + wy = cy + py⋅cy = cy(1 + py), and from these:
+;;
+;; lx⋅ly = cx⋅cy(1 - px)(1 - py)
+;; ux⋅uy = cx⋅cy(1 + px)(1 + py).
+;; ┌ ┐
+;; cx⋅cy(1 + px)(1 + py) - cx⋅cy(1 - px)(1 - py) │ factoring │
+;; Now, pz = ───────────────────────────────────────────── = │ out cx⋅cy │ =
+;; cx⋅cy(1 + px)(1 + py) + cx⋅cy(1 - px)(1 - py) └ ┘
+;; ┌ ┐
+;; (1 + px)(1 + py) - (1 - px)(1 - py) │ after expanding │ px + py
+;; = ─────────────────────────────────── = │ and │ = ─────────.
+;; (1 + px)(1 + py) + (1 - px)(1 - py) │ simplifying │ 1 + px⋅py
+;; └ ┘
+;; Assuming that px⋅py << 1, pz ≈ px + py. Test confirms our result:
+
+(define res1 (make-center-percent 4.7 0.02))
+(define res2 (make-center-percent 3.3 0.03))
+
+(define res1*res2 (mul-interval res1 res2))
+
+(percent res1) ; 0.02
+(percent res2) ; 0.03
+(percent res1*res2) ; 0.04997
View
52 2.1/2.14.scm
@@ -0,0 +1,52 @@
+
+;; Depends on some preceding exercises
+
+;; From the book
+(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)))))
+;; ----------------
+
+(define one (make-interval 1 1))
+
+;; A couple of resistors with 2% and 5% tolerance
+(define res3 (make-center-percent 39 0.02))
+(define res4 (make-center-percent 56 0.05))
+
+(center (par1 res3 res4)) ; 23.11
+(percent (par1 res3 res4)) ; 0.107
+
+(center (par2 res3 res4)) ; 22.98
+(percent (par2 res3 res4)) ; 0.032
+
+;; Results are different indeed, Lem is right
+
+;; Experiments with interval arithmetic
+(center (div-interval res3 res3)) ; 1.0008
+(percent (div-interval res3 res3)) ; 0.03998
+
+(center (div-interval res3 res4)) ; 0.6989
+(percent (div-interval res3 res4)) ; 0.0699 (dividing adds the uncertainties)
+
+(center (mul-interval res3 res4)) ; 2186 (should be 2184)
+(percent (mul-interval res3 res4)) ; 0.0699 (multiplying adds them too)
+
+(center (div-interval one res3)) ; 0.0257 (almost right)
+(percent (div-interval one res3)) ; 0.02
+(percent (div-interval one res4)) ; 0.05
+;; (taking reciprocal preserves the uncertainty)
+
+(center (add-interval res3 res4)) ; 95
+(percent (add-interval res3 res4)) ; 0.0377
+;; (uncertainty in adding is a bit more than the average of the tolerances)
+
+(center (sub-interval res4 res3)) ; 17
+(percent (sub-interval res4 res3)) ; 0.21
+;; (uncertainty in subtracting depends heavily on the difference --
+;; as difference gets smaller, relative uncertainty gets bigger)
View
47 2.1/2.15.scm
@@ -0,0 +1,47 @@
+
+;; Depends on some preceding exercises
+
+;; From the book
+(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)))))
+;; ----------------
+
+(define res3 (make-center-percent 39 0.05))
+(define res4 (make-center-percent 56 0.05))
+
+(center (par1 res3 res4)) ; 23.22
+(percent (par1 res3 res4)) ; 0.149
+
+(center (par2 res3 res4)) ; 22.99 (center is closer to actual value)
+(percent (par2 res3 res4)) ; 0.05 (par2 indeed makes tighter tolerances)
+
+(center (div-interval res3 res3)) ; 1.005 (should be exactly 1)
+(percent (div-interval res3 res3)) ; 0.0998 (should be 0)
+
+(center (div-interval res3 res4)) ; 0.6999
+(percent (div-interval res3 res4)) ; 0.0998 (dividing adds the uncertainties)
+
+;; She is right. We saw in 2.14 that different operations propagate relative
+;; tolerances differently - multiplication and division adds them,
+;; addition averages them and reciprocal preserves them. Subtraction usually
+;; scales them up.
+
+;; par1 has one multiplication, addition and division. We took the case where
+;; both resistances have the same percentage tolerance, 5%. Multiplication
+;; in numerator doubles the tolerance to 10%. Addition in denominator
+;; preserves tolerance at 5%. Dividing numerator by denominator adds the
+;; corresponding tolerances, producing 15%. Evaluating par1 confirms this:
+(percent (par1 res3 res4)) ; 0.149
+
+;; par2 has 1 addition and 3 reciprocals that all preserve the tolerance at 5%:
+(percent (par2 res3 res4)) ; 0.05
+
+;; Thus, we should prefer additions and reciprocals to multiplications and
+;; divisions.
View
134 2.1/2.16.scm
@@ -0,0 +1,134 @@
+
+;; The main reason why algebraically equivalent expressions behave differently
+;; in terms of uncertainty propagation is that the current system has "tunnel
+;; vision". It "sees" every arithmetic operation as separate from all others
+;; and performs them pairwise, producing new intervals. These are then again
+;; combined step-by-step with other intervals.
+
+;; The problem with this approach is that whenever there's a repeated
+;; occurrence of a variable in the formula, the program doesn't "remember"
+;; having seen it before. It treats this variable independently instead of
+;; realising that identically named variables deviate inside their intervals
+;; in lockstep, because they are obviously the "same thing".
+
+;; This explains why the formula with no repeated variables gives a better
+;; result with tighter interval. But it's not always possible to transform
+;; the formula so that every variable occurs just once. We need a different
+;; method that treats the formula as an atomic unit. One possibility is to
+;; find how sensitive the result is to the deviation of each variable,
+;; and then add every variable's contribution.
+
+;; "Small quantity" delta r
+(define Δr 0.1)
+
+;; Helper functions
+(define (nullify lst)
+ (map (λ (x) (* x 0)) lst))
+
+(nullify '(6 -3 0 12)) ; '(0 0 0 0)
+
+(define (shift-right lst)
+ (take (cons (last lst) lst) (length lst)))
+
+(define (shift-left lst)
+ (append (cdr lst) (list (car lst))))
+
+(shift-right '(1 0 0 0)) ; '(0 1 0 0)
+(shift-left '(1 0 0 0)) ; '(0 0 0 1)
+(shift-left '(0 0 0 1)) ; '(0 0 1 0)
+
+((repeated shift-right 3) '(1 0 0 0)) ; '(0 0 0 1)
+
+;; Combines arglist with Δr. Returns list of lists: first list is arglist
+;; where Δr is added to first element, second list is arglist where
+;; Δr is added to the second element, etc.
+(define (delcomb arglist)
+ (define startmask (cons Δr (nullify (rest arglist))))
+ (define (iter inlist outlist mask)
+ (if (empty? inlist)
+ outlist
+ (iter (rest inlist)
+ (append outlist (list (map + arglist mask)))
+ (shift-right mask))))
+ (iter arglist empty startmask))
+
+;; Delcomb with different implementation
+(define (delcomb arglist)
+ (define (iter inlist auxlist outlist)
+ (if (empty? inlist)
+ outlist
+ (iter (rest inlist)
+ (append auxlist (list (first inlist)))
+ (append outlist
+ (list (append auxlist
+ (list (+ (first inlist) Δr))
+ (rest inlist)))))))
+ (iter arglist empty empty))
+
+;; This takes in quoted function f and list of intervals (arglist).
+;; Arity of f and number of items in arglist must be the same.
+;; Returns an interval whose center is calculated by applying f to the
+;; centers of argument intervals. Uncertainty of the returned value
+;; is found by calculating how much each argument's uncertainty contributes
+;; to the result's uncertainty, and adding these contributions.
+(define (f-interval f arglist)
+ (let* ((fargs (second f))
+ (centers (map center arglist))
+ (widths (map width arglist))
+ (centers+Δr (delcomb centers)))
+ (if (not (= (length fargs) (length arglist)))
+ (printf "Arity of the function and arglist length must match.\n")
+ (let* ((f-center (eval-f f centers))
+ (f-width (apply + (map (λ (center+Δr width)
+ (* (/ (- (eval-f f center+Δr)
+ f-center)
+ Δr)
+ width))
+ centers+Δr widths))))
+ (make-center-width f-center f-width)))))
+
+;; Treats interval widths as standard deviations and adds them in quadrature
+(define (f-interval f arglist)
+ (let* ((fargs (second f))
+ (centers (map center arglist))
+ (widths (map width arglist))
+ (centers+Δr (delcomb centers)))
+ (if (not (= (length fargs) (length arglist)))
+ (printf "Arity of the function and arglist length must match.\n")
+ (let* ((f-center (eval-f f centers))
+ (f-width (sqrt
+ (apply + (map (λ (center+Δr width)
+ (sqr (* (/ (- (eval-f f center+Δr)
+ f-center)
+ Δr)
+ width)))
+ centers+Δr widths)))))
+ (make-center-width f-center f-width)))))
+
+(define (eval-f f arglist)
+ (eval (append (list f) arglist)))
+
+;; Parallel resistance formula
+(define par3 '(λ (R1 R2) (/ (* R1 R2) (+ R1 R2))))
+
+;; Equivalent to previous formula
+(define par4 '(λ (R1 R2) (/ (+ (/ R1) (/ R2)))))
+
+;; Three parallel resistances
+(define par5 '(λ (R1 R2 R3) (/ (+ (/ R1) (/ R2) (/ R3)))))
+
+;; We take two resistors...
+(define res3 (make-center-percent 39 0.05))
+(define res4 (make-center-percent 56 0.05))
+
+;; ...and connect them in parallel...
+(center (f-interval par3 (list res3 res4))) ; 22.989
+(percent (f-interval par3 (list res3 res4))) ; 0.0499
+
+;; ...and again using equivalent formula (should give same result)
+(center (f-interval par4 (list res3 res4))) ; 22.989
+(percent (f-interval par4 (list res3 res4))) ; 0.0499 (same indeed)
+
+;; Should be third of res3
+(center (f-interval par5 (list res3 res3 res3))) ; 13.0
+(percent (f-interval par5 (list res3 res3 res3))) ; 0.0499

0 comments on commit cc66732

Please sign in to comment.