Dan's Thoughts Thinking somewhat carefully

19May/100

sin(n) = -1 where n is an integer in radians

I saw a fun math problem on reddit the other day.

find a number n, so that sin(n)=-1, where n is an integer in radians; so sin(270 degrees) doesn't count. Obviously it will never be exactly -1 but close enough for the difference to be lost in rounding.

\sin(\frac{3\pi}{2} + 2\pi \cdot n) = -1
I need the argument of sin function to be as close to an integer as possible. Call this integer m.
\frac{3\pi}{2}+2\pi \cdot n \approx m
Solving for \pi leads to:
\pi \approx \frac{p}{q} \approx \frac{2m}{4n+3}

If I have a rational approximation to \pi with an even numerator I can divide it by two get my m. I also have to make sure that the denominator is in the form of 4n+3.
It's possible to use continued fractions to approximate real numbers. Here's a continued fraction sequence for pi: http://www.research.att.com/~njas/sequences/A001203

The first rational approximation I learned in elementary school is 22/7 which is perfect.

> (sin 11)
-.9999902065507035

For the others I'll have to evaluate the continued fraction to get my approximation of a simple fraction.

> (eval-terms (list 3 7 15 1 292 1))
104348/33215

> (sin (/ 104348 2))
-.9999999999848337

> (eval-terms (list 3 7 15 1 292 1 1 1 2 1 3 1 14 2 1))
245850922/78256779

> (sin (/ 245850922 2))
-.99999999999999999532
Looks like a good candidate was found.

This is the code to evaluate a continued fraction coefficients. It's very convenient that scheme has a native rational data type.

(define (eval-terms ts)
  (cond ((= (length ts) 1) (car ts))
        (else
          (+ (car ts) (/ 1 (eval-terms (cdr ts)))))))
 
(define (eval-terms-iter ts)
  (let ((rev-ts (reverse ts)))
    (define (ev-terms ts frac)
      (if (null? ts)
        frac
        (ev-terms (cdr ts) (+ (car ts) (/ 1 frac)))))
    (ev-terms rev-ts (car rev-ts))))
Filed under: lisp, math No Comments
21Oct/090

QMISMF Chap 3

3-1
\begin{pmatrix} 1 & 0 \\ 0 & -1  \end{pmatrix}+\begin{pmatrix} 0 & -i \\ i & 0  \end{pmatrix}=\begin{pmatrix} 1 & -i \\ i & -1  \end{pmatrix}
\begin{pmatrix} 0 & 1 \\ 1 & 0  \end{pmatrix}+\begin{pmatrix} 0 & -i \\ i & 0  \end{pmatrix}=\begin{pmatrix} 0 & 1-i \\ 1+i & 0 \end{pmatrix}
\frac{1}{2}\begin{pmatrix} 1 & 0 \\ 0 & 1  \end{pmatrix}+\frac{1}{2}\begin{pmatrix} 1 & 0 \\ 0 & -1  \end{pmatrix}=\frac{1}{2}\begin{pmatrix} 2 & 0 \\ 0 & 0 \end{pmatrix}=\begin{pmatrix} 1 & 0 \\ 0 & 1  \end{pmatrix}
\begin{pmatrix} 0 & 1 \\ 1 & 0  \end{pmatrix}+\begin{pmatrix} 0 & 1 \\ -1 & 0  \end{pmatrix}=\begin{pmatrix} 0 & 2 \\ 0 & 0  \end{pmatrix}
3-2
\begin{pmatrix} 1 & -1 & 1 \\ 0 & 1 & 0 \\ 2 & 0 & 3  \end{pmatrix} \begin{pmatrix} 3 & 3 & -1 \\ 0 & 1 & 0 \\ -2 & -2 & 1  \end{pmatrix}=\begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1  \end{pmatrix}
\begin{pmatrix} 3 & 3 & -1 \\ 0 & 1 & 0 \\ -2 & -2 & 1  \end{pmatrix} \begin{pmatrix} 1 & -1 & 1 \\ 0 & 1 & 0 \\ 2 & 0 & 3  \end{pmatrix} =\begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1  \end{pmatrix}
3-3
\frac{1}{4}\begin{pmatrix} 1 & 1 \\ 1 & 1  \end{pmatrix}\begin{pmatrix} 1 & 1 \\ 1 & 1  \end{pmatrix}=\frac{1}{4}\begin{pmatrix} 2 & 2 \\ 2 & 2  \end{pmatrix}=\frac{1}{2}\begin{pmatrix} 1 & 1 \\ 1 & 1  \end{pmatrix}
\frac{1}{2}\begin{pmatrix} 1 & -1 \\ -1 & 1  \end{pmatrix}\frac{1}{2}\begin{pmatrix} 1 & -1 \\ -1 & 1  \end{pmatrix}=\frac{1}{4}\begin{pmatrix} 2 & -2 \\ -2 & 2  \end{pmatrix}=\frac{1}{2}\begin{pmatrix} 1 & -1 \\ -1 & 1  \end{pmatrix}
\frac{1}{2}\begin{pmatrix} 1 & 1 \\ 1 & 1  \end{pmatrix}\frac{1}{2}\begin{pmatrix} 1 & -1 \\ -1 & 1  \end{pmatrix}=\frac{1}{4}\begin{pmatrix} 0 & 0 \\ 0 & 0  \end{pmatrix}=\begin{pmatrix} 0 & 0 \\ 0 & 0  \end{pmatrix}
\frac{1}{2}\begin{pmatrix} 1 & -1 \\ -1 & 1  \end{pmatrix}\frac{1}{2}\begin{pmatrix} 1 & 1 \\ 1 & 1  \end{pmatrix}=\frac{1}{4}\begin{pmatrix} 0 & 0 \\ 0 & 0  \end{pmatrix}=\begin{pmatrix} 0 & 0 \\ 0 & 0  \end{pmatrix}
3-4
\begin{pmatrix} 1 & 2 & 3 \\ 2 & 4 & 5 \\ 3 & 5 & 6  \end{pmatrix} \begin{pmatrix} 1 & -3 & 2 \\ -3 & 3 & -1 \\ 2 & -1 & 0 \end{pmatrix} =\begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}
3-5
Assume: AM^{-1}=M^{-1}A \iff AM=MA
AM^{-1}=M^{-1}A \implies AM=MA
AM^{-1}=M^{-1}A
AM^{-1}M=M^{-1}AM
MAI=MM^{-1}AM
MA=IAM
MA=AM

AM=MA \implies AM^{-1}=M^{-1}A
AM=MA
M^{-1}AM=A
M^{-1}A=AM^{-1}
3-6
Assume: AM=MA and BM=MB
(A+B)M=M(A+B)
AM+BM=MA+MB
AM+BM-MA-MB=0
(AM-MA)+(BM-MB)=0
0+0=0 (AB)M=M(AB)
A(BM)=(MA)B
A(MB)=(AM)B
(AM)B=(AM)B (zB)M=M(zB)
z(BM)=(MB)z
z(BM)=(BM)z

3-7
\begin{pmatrix} 2 & 2-i  \\ 2+i & -2 \end{pmatrix} \frac{1}{9}\begin{pmatrix} 2 & 2-i  \\ 2+i & -2 \end{pmatrix}=\frac{1}{9}\begin{pmatrix} 9 & 0  \\ 0 & 9 \end{pmatrix}=\begin{pmatrix} 1 & 0  \\ 0 & 1 \end{pmatrix}

Filed under: math, physics, qmismf No Comments
5Oct/090

QMISMF: Chapter 2

2-1:

> (+ 3-i 2+4i)
5+3i
> (+ 1+3i 2)
3+3i
> (- -5+2i 2+2i)
-7
> (+ -2+i 2+2i)
+3i
> (* 3-i 2+4i)
10+10i
> (* 1+3i 2)
2+6i
> (* 0+i 1+3i)
-3+i
> (* -5+2i 2+3i)
-16-11i
> (* 2+3i -2+3i)
-13
> (* 2+3i 3+2i)
+13i

2-2:

z^{-1}=\frac{1}{x^2+y^2}(x+iy)
\left(\frac{1}{i}\right)^{-1} \rightarrow 0-i=\frac{1}{1}(i)=i
i \cdot \frac{1}{i}= 1
\left(\frac{1}{-i}\right)^{-1} \rightarrow \frac{-i}{-1} = i = \frac{1}{1}(-i)=-i
 -i \cdot \frac{1}{-i}= 1

2-3:

i^*i = -i \cdot i = 1
(-i)^*(-i)= i \cdot -i = 1
|i|^2=(i^*i)=1
|-i|^2=\left( (-i)^*(-i) \right)=1
|i| = \sqrt{i^*i} = \sqrt{-i \cdot i} = 1
|-i| = \sqrt{(-i)^*(-i)}= \sqrt{i \cdot -i} = 1

2-4:

z=3+4i
z^* = 3-4i
z^*z = (3-4i)(3+4i) = 25
|z| = \sqrt{z^*z} = \sqrt{25} = 5
z^2 = z \cdot z = (3+4i)(3+4i) = -7+24i
\frac{1}{z} = \frac{1}{z} \cdot \frac{z^*}{z^*} = \frac{z^*}{zz^*} = \frac{3-4i}{25}

2-5:

|w+z| < |w| + |z|
|3+4i| < |3| + |4i|
5 < 7

2-6:

\sqrt{z} = \sqrt{\frac{\sqrt{x^2+y^2}}{2}} \cdot \left( \sqrt{1 + \frac{x}{\sqrt{x^2+y^2}}} + i\sqrt{1 - \frac{x}{\sqrt{x^2+y^2}}} \right)
for y \ge 0
\sqrt{z} = \sqrt{\frac{\sqrt{x^2+y^2}}{2}} \cdot \left( -\sqrt{1 + \frac{x}{\sqrt{x^2+y^2}}} + i\sqrt{1 - \frac{x}{\sqrt{x^2+y^2}}} \right)
for y \le 0 \sqrt{i} = \sqrt{0+1i} = \sqrt{\frac{\sqrt{0^2+1^2}}{2}} \cdot \left( \sqrt{1 + \frac{0}{\sqrt{0^2+1^2}}} + i\sqrt{1 - \frac{0}{\sqrt{0^2+1^2}}} \right)
\sqrt{i} = \sqrt{0+1i} = \sqrt{\frac{1}{2}} \cdot (1+i) = \sqrt{\frac{1}{2}}+i\sqrt{\frac{1}{2}}
\left(\sqrt{\frac{1}{2}}+i\sqrt{\frac{1}{2}}\right)^2 = \frac{1}{2} + i - \frac{1}{2} = i \sqrt{-i} = \sqrt{0-1i} = \sqrt{\frac{\sqrt{0^2+(-1)^2}}{2}} \cdot \left( -\sqrt{1 + \frac{0}{\sqrt{0^2+(-1)^2}}} + i\sqrt{1 - \frac{0}{\sqrt{0^2+(-1)^2}}} \right)
\sqrt{-i} = \sqrt{0-1i} = \sqrt{\frac{1}{2}} \cdot (-1+i) = -\sqrt{\frac{1}{2}}+i\sqrt{\frac{1}{2}}
\left(- \sqrt{\frac{1}{2}}+i \sqrt{\frac{1}{2}}\right) \cdot \left(- \sqrt{\frac{1}{2}}+i \sqrt{\frac{1}{2}}\right)=\frac{1}{2}-i-\frac{1}{2}=-i

2-7:

z=-\frac{b}{2a} \pm \sqrt{\left(\frac{b}{2a}\right)^2 - \frac{c}{a}}
z=-\frac{b}{2a} \pm \sqrt{\frac{b^2-4ac}{2^2a^2}}
z=-\frac{b}{2a} \pm \frac{1}{2a}\sqrt{b^2-4ac}
z=\frac{-b \pm \sqrt{b^2-4ac}}{2a} since the above is now in the usual form of the quadratic formula it is clear z is a solution.

Filed under: lisp, math, physics, qmismf No Comments
23Aug/090

Exercise 2.38 of SICP

Exercise 2.38: The accumulate procedure is also known as fold-right, because it combines the first element of the sequence with the result of combining all the elements to the right. There is also a fold-left, which is similar to fold-right, except that it combines elements working in the opposite direction:

(define (fold-left op initial sequence)
  (define (iter result rest)
    (if (null? rest)
        result
        (iter (op result (car rest))
              (cdr rest))))
  (iter initial sequence))

What are the values of

(fold-right / 1 (list 1 2 3))

3/2

(fold-left / 1 (list 1 2 3))

1/6

(fold-right list nil (list 1 2 3))

(1 (2 (3 ())))

(fold-left list nil (list 1 2 3))

(((() 1) 2) 3)
Give a property that op should satisfy to guarantee that fold-right and fold-left will produce the same values for any sequence.
(fold-right op i (list a1 a2 a3))
(a_1 \cdot (a_2 \cdot ( a_3 \cdot i)))
(fold-left op i (list a1 a2 a3))
(((i \cdot a_1) \cdot a_2) \cdot a_3)

Any binary associative operation will be invariant under fold-right and fold-left.

Since 2.37 involved matrix multiplication which is associative, I will use that as an example.
> (define i (list (list 1 0 0) (list 0 1 0) (list 0 0 1)))
> (define a1 (list (list 8 3 2) (list 1 0 9) (list 3 4 5)))
> (define a2 (list (list 5 6 7) (list 1 2 8) (list 6 7 7)))
> (define a3 (list (list 6 7 9) (list 4 3 1) (list 3 4 5)))
> (fold-left matrix-*-matrix i (list a1 a2 a3))
((884 965 1033) (840 900 950) (802 878 942))
> (fold-right matrix-*-matrix i (list a1 a2 a3))
((884 965 1033) (840 900 950) (802 878 942))

Notice that i is the identity matrix because if it wasn't, I would also be testing commutativity of matrix product. Matrix product doesn't commute. If i is changed to another matrix fold-left and fold-right will return differing results.

Filed under: SICP, lisp, math No Comments
22Aug/090

Exercise 2.37 of SICP

Exercise 2.37: Suppose we represent vectors v=(v_i) as sequences of numbers, and matrices m = (m_{ij}) as sequences of vectors (the rows of the matrix). For example, the matrix

M = \begin{bmatrix} 1&2&3&4\\ 4&5&6&6\\ 6&7&8&9 \end{bmatrix}

is represented as the sequence ((1 2 3 4) (4 5 6 6) (6 7 8 9)). With this representation, we can use sequence operations to concisely express the basic matrix and vector operations. These operations (which are described in any book on matrix algebra) are the following:

(dot-product v w)
returns the sum \sum_{i} v_iw_i
(matrix-*-vector m v)
returns the vector t, where t_i = \sum_{j} m_{ij}v_j
(matrix-*-matrix m n)
returns the matrix p where, p_{ij} = \sum_{k} m_{ik}n_{kj}
(transpose m)
returns the matrix n where, n_{ij} = m_{ji}

We can define the dot product as:

(define (dot-product v w)
  (accumulate + 0 (map * v w)))

Fill in the missing expressions in the following procedures for computing the other matrix operations. (The procedure accumulate-n is defined in exercise 2.36.)

(define (matrix-*-vector m v)
  (map <??> m))
(define (transpose mat)
  (accumulate-n <??> <??> mat))
(define (matrix-*-matrix m n)
  (let ((cols (transpose n)))
    (map <??> m)))

Solutions:

(define (accumulate fn init-value items)
  (if (null? items)
    init-value
    (fn (car items)
        (accumulate fn init-value (cdr items)))))
 
(define (accumulate-n op init seqs)
  (if (null? (car seqs))
    '()
      (cons (accumulate op init (map (lambda (e) (car e)) seqs))
            (accumulate-n op init (map (lambda (e) (cdr e)) seqs)))))
(define (matrix-*-vector m v)
  (map (lambda (row) (dot-product row v)) m))
(define (transpose mat)
  (accumulate-n cons '() mat))
(define (matrix-*-matrix m n)
  (let ((cols (transpose n)))
    (map (lambda (row) (matrix-*-vector cols row))
         m)))

> (define n (list (list 1 2) (list 3 4)))
> (matrix-*-vector n (list 10 20))
(50 110)
> (transpose n)
((1 3) (2 4))
> (define m (list (list 1 2 3) (list 4 5 6) (list 7 8 9)))
> (matrix-*-matrix m m)
((30 36 42) (66 81 96) (102 126 150))
> (define q (list (list 9 1) (list 3 4) (list 6 7) (list 1 8)))
> (matrix-*-matrix (list (list 1 2 3 4) (list 5 6 7 8)) q)
((37 62) (113 142))

Filed under: SICP, lisp, math No Comments
19Aug/090

Exercise 2.34 of SICP

Exercise 2.34: Evaluating a polynomial in x at a given value of x can be formulated as an accumulation. We evaluate the polynomial
p(x) = a_n x^n + a_{n-1}x^{n-1} + \cdots + a_1 x+ a_0
using a well-known algorithm called Horner's rule, which structures the computation as
(\cdots (a_n x + a_{n-1})x + \cdots + a_1)x + a_0
In other words, we start with an, multiply by x, add an-1, multiply by x, and so on, until we reach a0. Fill in the following template to produce a procedure that evaluates a polynomial using Horner's rule. Assume that the coefficients of the polynomial are arranged in a sequence, from a0 through an.

(define (horner-eval x coefficient-sequence)
  (accumulate (lambda (this-coeff higher-terms) <??>)
              0
              coefficient-sequence))

For example, to compute 1 + 3x + 5x^3 + x^5 at x = 2 you would evaluate

(horner-eval 2 (list 1 3 0 5 0 1))

I included a non-abstracted version of Horner's method for clarity.

(define (accumulate fn init-value items)
  (if (null? items)
    init-value
    (fn (car items)
        (accumulate fn init-value (cdr items)))))
 
(define (he x coeff-seq)
  (if (null? coeff-seq)
    0
    (+ (car coeff-seq)
       (* x (he x (cdr coeff-seq))))))
 
(define (horner-eval x coeff-seq)
  (accumulate (lambda (this-coeff higher-terms)
                (+ this-coeff (* x higher-terms)))
              0
              coeff-seq))

> (horner-eval 2 (list 1 3 0 5 0 1))
79
> (horner-eval 2 (list 1))
1
> (horner-eval 2 (list ))
0
> (horner-eval 2 (list 1 1))
3

Filed under: SICP, lisp, math No Comments
17Aug/090

Exercise 2.32 of SICP

Exercise 2.32: We can represent a set as a list of distinct elements, and we can represent the set of all subsets of the set as a list of lists. For example, if the set is (1 2 3), then the set of all subsets is (() (3) (2) (2 3) (1) (1 3) (1 2) (1 2 3)). Complete the following definition of a procedure that generates the set of subsets of a set and give a clear explanation of why it works:

(define (subsets s)
  (if (null? s)
      (list nil)
      (let ((rest (subsets (cdr s))))
        (append rest (map <??> rest)))))

In order for me to understand what is going on here, I understand how powersets are constructed first and then derive this function.

This program wants to create what is called a powerset.
If S is a set of {a1,a2,a3}
\mathcal{P}(S) = \{\{\} \{a_3\} \{a_2\} \{a_2 a_3\} \{a_1\} \{a_1 a_3\} \{a_1 a_2\} \{a_1 a_2 a_3\}\}
The reasoning to construct a powerset is similar to the way the count change example is thought about.
Create a powerset for the cdr elements without using the car element, cons the car element with powerset of cdr elements. To the result append powerset of cdr elements.

When the arguments are an element and a set inject the element into every element of the set creating a new set. Append that result with the original powerset injected into.

P{\varnothing}
{\varnothing}

P{a1}
(append (inject a1 P{\varnothing}) P{\varnothing})
{{a1} \varnothing}

P{a1 a2}
(append (inject a1 P{a2}) P{a2})
(append {{a1 a2} {a1}} P{a2})
{{a1 a2} {a1} {a2} \varnothing}

P{a1 a2 a3}
(append (inject a1 P{a2 a3}) P{a2 a3})
(append (inject a1 (apend (inject a2 P{a3}) P{3}))
        (append (inject a2 P{a3}) P{a3}))

(append (inject a1 (append (inject a2 {{a3} \varnothing}) {{a3} \varnothing}))
        (append (inject a2 {{a3} \varnothing}) {{a3} \varnothing}))

(append (inject a1 (append {{a3 a2} {a2}} {{a3} \varnothing}))
        (append (inject a2 {{a3} \varnothing}) {{a3} \varnothing}))

(append (inject a1 {{a3 a2} {a2} {a3} \varnothing})
        (append {{a3 a2} {a2}} {{a3} \varnothing}))

(append {{a3 a2 a1 } {a2 a1} {a3 a1} {a1}}
        (append {{a3 a2} {a2}} {{a3} \varnothing}))

(append {{a3 a2 a1 } {a2 a1} {a3 a1} {a1}}
        {{a3 a2} {a2} {a3} \varnothing})

{{a3 a2 a1} {a2 a1} {a3 a1} {a1} {a3 a2} {a2} {a3} \varnothing}

Here is the process above in code:

(define (powerset ss)
  (if (null? ss)
    (list '())
    (append (inject (car ss) 
                    (powerset (cdr ss))) 
            (powerset (cdr ss)))))
(define (inject s ss)
  (map (lambda (elem)
           (cons s elem))
       ss))

Finally:

(define (subsets s)
  (if (null? s)
      (list '())
      (let ((rest (subsets (cdr s))))
        (append rest (map (lambda (elem) (cons (car s) elem)) 
                          rest)))))

It's worth noticing that because of the let expression (subsets (cdr s)) is only called once, as opposed to twice in my earlier derivation.

Filed under: SICP, lisp, math No Comments
1Aug/090

Exercise 2.16 of SICP

Exercise 2.16: Explain, in general, why equivalent algebraic expressions may lead to different answers. Can you devise an interval-arithmetic package that does not have this shortcoming, or is this task impossible? (Warning: This problem is very difficult.)

The errors are a lot tighter for par2 than par1. The reason that this is true is due to the fact that any operation between two intervals will increase errors. If two algebraically equivalent statements are evaluated, the statement with the least operations between intervals will produce less errors.

Here's a slightly exaggerated example:
A = 10 \pm 0.1
A = \frac{A^7}{A^6} = \frac{A \cdot A \cdot A \cdot A \cdot A \cdot A \cdot A}{A \cdot A \cdot A \cdot A \cdot A \cdot A}

> (define A (make-center-percent 10 1))
> (print A)
[10.,.9999999999999963]
> (print (div-interval (interval-pow A 7) (interval-pow A 6)))
[10.084120477031101,12.92768447766068]

The most noticeable difference is between percent errors of 1% and 13%. Even though the fraction is algebraically equivalent to A evaluating the intervals gives very different answers.

(define (interval-pow x n)
  (define (square i)
    (mul-interval i i))
  (cond ((= n 1) x)
        ((= 0 (remainder n 2)) (square (interval-pow x (/ n 2))))
        ((= 1 (remainder n 2)) (mul-interval x (interval-pow x (- n 1))))))

In order to devise a system that doesn't have this problem would mean it would have to be smart enough to detect duplicate entries and rewrite them in a more interval arithmetic efficient manner. At a minimum it would have to analyze code of this form:

(div-interval (interval-pow A 7) (interval-pow A 6))

and reduce it to A and than evaluate it.

Transforming the resistor example is even tougher.

(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)))))

The program would have to recognize division and simplify it by dividing the numerator and denominator by r1 and r2. Both the exponent and parallel resistor examples would have to be transformed symbolically since simplifying arithmetically would only introduce errors. I do not think it is feasible to write a program like this with the current lisp knowledge learned in the book so far.

Filed under: SICP, lisp, math No Comments
31Jul/090

Exercise 2.15 of SICP

Exercise 2.15: Eva Lu Ator, another user, has also noticed the different intervals computed by different but algebraically equivalent expressions. She says that a formula to compute with intervals using Alyssa's system will produce tighter error bounds if it can be written in such a form that no variable that represents an uncertain number is repeated. Thus, she says, par2 is a "better'' program for parallel resistances than par1. Is she right? Why?

> (print (div-interval (make-center-percent 6.8 1) (make-center-percent 6.8 1)))
[1.0002000200020003,1.9998000199979908]
> (print (div-interval (make-center-percent 6.8 .1) (make-center-percent 6.8 .1)))
[1.000002000002,.19999980000021655]

> (print (div-interval (make-center-percent 6.8 .1) (make-center-percent 3.4 .1)))
[2.000004000004,.19999980000021655]
> (print (div-interval (make-center-percent 6.8 .1) (make-center-percent 1.7 .1)))
[4.000008000008,.19999980000021655]
> (print (div-interval (make-center-percent 6.8 .1) (make-center-percent 1.7 .01)))
[4.000000440000004,.1099999890000035]

> (print (par1 (make-center-percent 6.8 .1) (make-center-percent 1.7 .01)))
[1.3600022771855311,.19199981581619266]
> (print (par2 (make-center-percent 6.8 .1) (make-center-percent 1.7 .01)))
[1.3599998237438813,.028000014256019334]

> (print (par1 (make-center-percent 6.8 .1) (make-center-percent 6.8 .1)))
[3.4000136000136,.2999992000024115]
> (print (par2 (make-center-percent 6.8 .1) (make-center-percent 6.8 .1)))
[3.4000000000000004,.10000000000000857]

> (print (par1 (make-center-percent 6.8 10) (make-center-percent 4.7 5)))
[2.844199964577264,22.613352145193346]
> (print (par2 (make-center-percent 6.8 10) (make-center-percent 4.7 5)))
[2.777440701636504,7.05260392723452]

(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)))))

Alyssa's system will produce tighter error bounds. The errors are a lot tighter for par2 than par1. The reason that this is true is due to the fact that any operation between two intervals will increase errors. If two algebraically equivalent statements are evaluated, the statement with the least operations between intervals will produce less errors.

Here's a slightly exaggerated example:
A = 10 \pm 0.1
A = \frac{A^7}{A^6} = \frac{A \cdot A \cdot A \cdot A \cdot A \cdot A \cdot A}{A \cdot A \cdot A \cdot A \cdot A \cdot A}

> (define A (make-center-percent 10 1))
> (print A)
[10.,.9999999999999963]
> (print (div-interval (interval-pow A 7) (interval-pow A 6)))
[10.084120477031101,12.92768447766068]

The most noticeable difference is between percent errors of 1% and 13%.

(define (interval-pow x n)
  (define (square i)
    (mul-interval i i))
  (cond ((= n 1) x)
        ((= 0 (remainder n 2)) (square (interval-pow x (/ n 2))))
        ((= 1 (remainder n 2)) (mul-interval x (interval-pow x (- n 1))))))
Filed under: SICP, lisp, math No Comments
29Jul/090

Exercise 2.13 of SICP

Exercise 2.13: Show that under the assumption of small percentage tolerances there is a simple formula for the approximate percentage tolerance of the product of two intervals in terms of the tolerances of the factors. You may simplify the problem by assuming that all numbers are positive.

(a-\epsilon_1,a+\epsilon_1) \cdot (b-\epsilon_2,b+\epsilon_2)
(a+\epsilon_1)\cdot(b+\epsilon_2)=ab+a\epsilon_2+b\epsilon_1+\epsilon_1\epsilon_2
(a+\epsilon_1)\cdot(b-\epsilon_2)=ab-a\epsilon_2+b\epsilon_1-\epsilon_1\epsilon_2
(a-\epsilon_1)\cdot(b+\epsilon_2)=ab+a\epsilon_2-b\epsilon_1-\epsilon_1\epsilon_2
(a-\epsilon_1)\cdot(b-\epsilon_2)=ab-a\epsilon_2-b\epsilon_1+\epsilon_1\epsilon_2

Since a>0 and b>0:
(a+\epsilon_1)\cdot(b+\epsilon_2) \approx ab+ a\epsilon_2+b\epsilon_1+\epsilon_1\epsilon_2
(a-\epsilon_1)\cdot(b-\epsilon_2) \approx ab -a\epsilon_2-b\epsilon_1+\epsilon_1\epsilon_2
Assuming \epsilon_1 and \epsilon_2 are small, the quantity \epsilon_1\epsilon_2 can be discarded.
(a+\epsilon_1)\cdot(b+\epsilon_2) \approx ab+a\epsilon_2+b\epsilon_1
(a-\epsilon_1)\cdot(b-\epsilon_2) \approx ab-a\epsilon_2-b\epsilon_1

Since we are interested in the tolerance terms:
(a+\epsilon_1)\cdot(b+\epsilon_2) \approx a\epsilon_2+b\epsilon_1
(a-\epsilon_1)\cdot(b-\epsilon_2) \approx -a\epsilon_2-b\epsilon_1

The new product interval in terms of small tolerance terms looks like this:
[ab-(a\epsilon_2+b\epsilon_1),ab+(a\epsilon_2+b\epsilon_1)]

Filed under: SICP, math No Comments