You are currently browsing the tag archive for the ‘numerical analysis’ tag.

Procedures are like mathematical functions, with one important difference: a mathematical function can provide declarative definitions, whereas a procedure must provide imperative definitions. A declarative statement is a ‘this is true’ statement. An imperative statement is a ‘do this’ statement. We’ll examine the difference by looking at how the square root function is defined in mathematics and in computer science.

### 1.1.7 Example: Square roots by Newton’s method

When we write sqrt(x) in mathematics, what do we mean? One possible definition is

sqrt(x) is the unique non-negative real number whose square is equal to x

or, expressed in mathematical notation,

$\sqrt{x}$ is the unique $y\in\mathbb{R}$ such that $y\geq 0$ and $y^2=x$.

This is a perfectly valid definition, and you can prove that it is well-formed: that is, you can prove that for non-negative numbers x, there really is a unique non-negative number y such that the square of y is x, and therefore that the definition of the ‘sqrt’ function makes sense.

However, just knowing that we have a sensibly defined function is not much help in computer science. We need to know how to computer the function. That is, we need a sequence of statements that that lets us deduce the value of the output of the function, given an input (actually we will only compute an approximation to the value of the function, but that’s okay – we can make our approximations really very accurate).

This is what we mean by declarative vs imperative knowledge. According to SICP, mathematics is concerned primarily with declarative (what is) statements, whereas computer science is primarily concered with imperative (how to) statements. I’d take issue with this, given the existence of a whole subfield of mathematics, numerical analysis, which concerns itself with how best to find solutions to equations. But I’m not here to debate semantics – I’m here to learn about computer science.

So how do we compute square roots?

The most common way is to use Newton’s method of successive approximations (tip: if you’re interviewing for a job that involves mathematical and programming literacy, learn how to compute square roots using Newton’s method. It’s a favourite interview question). Newton’s method is a formula for finding roots of general functions $f$, but here we only need it in a very simple form. To find the solution of y = sqrt(x), we make an initial guess (say y = 1) and check if it’s accurate enough. If not, then we refine our guess by averaging y and x/y, and check again. We continue iterating in this way until we eventually converge on the solution. You can prove that for any positive initial guess, this method converges to the right solution. Moreover, you can prove that you double the number of digits of accuracy with each iteration — this is known as quadratic convergence.

Let’s formalize this in terms of procedures. This is just a matter of translating our description of the formula above into Lisp code:

(define (sqrt-iter guess x)
(if (good-enough? guess x)
guess
(sqrt-iter (improve guess x)
x)))


This procedure simply says that, given a guess for the solution, we first check if the guess is good enough; if it is, then we report the guess as our solution, if not, then we improve the guess and apply the procedure again. Note that we haven’t yet defined the procedures good-enough? or improve yet. We’re using a programming strategy called wishful thinking — we write the program as if the sub-procedures we need exist, and then we write them later.

We already know how to improve our guess: we average the old guess y with x/y, like this:

(define (improve guess x)
(average guess (/ x guess)))


where the procedure average is defined to be

(define (average x y)
(/ (+ x y) 2))


We also have to define what we mean by ‘good enough’. We’ll say that our guess is good enough if the square of the guess is within a small distance (say 0.001) of the input. This isn’t a very good test, for several reasons, but it will do for now:

(define (good-enough? guess x)
(< (abs (- (square guess) x)) 0.001))


Finally, we need a starting guess. We can always guess that the square root of a number is 1:

(define (sqrt x)
(sqrt-iter 1.0 x))


Note that because Scheme has a rational data type, we need to put 1.0 as our initial guess, rather than 1. If x was an integer and we guessed 1, then the interpreter would use rational arithmetic, and all subsequence operations would result in rational numbers. By guessing 1.0 instead, we force the result of subsequent operations to be floating point numbers.

We can now use our new procedure as we’d use any other:

> (sqrt 9)
3.00009155413138

> (sqrt 100)
10.000000000139897

> (square (sqrt 1000))
1000.000369924366


This demonstrates that the simple procedural language we have introduced so far is suffient for writing numerical programs. Notice that we haven’t introduced any iterative (looping) constructes yet. Instead, we relied on the simple ability of a procedure to call itself.