Tuesday, February 12, 2008
I'm a teacher's assistant for the Physics Department's Numerical Methods course. The course is taught in Matlab (with some option "bonus" material handling C or Fortran) and means to take the student from root finding to partial differential equations, visiting subjects like polynomial approximation in between. It is that latter subject which I wish to discuss today, particularly the notion of divided differences. Suppose you wish to approximate a function's form after having been given a set of pairs (xi,fi). It turns out that the following quantities are very useful in writing the coefficients of a polynomial (whose degree is commensurate with the number of points you have). DD1,2(x,f) = f2-f1/(x2-x1) with subsequent differences as DDi...j(x,f) = DDi+1...j(x,f)-DDi...j-1(x,f)/(xj-xi) Most readers with a bit of experience programming will immediately, upon being asked to calculate these numbers, lock upon the obvious recursive nature of the problem and might implement a solution like this (in Scheme - since I don't want my students lifting the code from my blog).
Which seems simple enough (if you are have bothered to learn scheme at this point and can peer through the parenthesis). Notably, recursive definition is almost a transcription of the definition (which is in the book). Implementing the calculation of the divided difference directly, with a set of while and for loops, is comparatively more complex. Nevertheless, after some argument, we agreed not to attempt to instruct the students in the style of recursion. I felt that anyone who wants to claim to be familiar with programming, even somewhat superficially, should understand recursion enough to translate obviously recursive concepts into recursive code. Its true that the above code is not tail-recursive (writing a tail recursive version necessitates an increase in complexity which is beyond the scope of this entry) but for the typical use case (a few points, no more than 20 at the absolute most), the recursive solution is clear, concise and simple to debug and it offers acceptable performance for a course which has sworn off efficiency as a pedagogical goal. The reason, according to the Professor in charge, is that people do not sit down with a pen and paper to perform a calculation and proceed "recursively." And it is true, if you asked me to calculate a particular divided difference I would start at the bottom and work my way up, serially (aware of what lower differences I needed to compute higher ones). I think this accounts for the "backwards" or "inside out" feeling people have towards recursion: we think and operate with a global scope when we sit down to calculate (say, the page we work on) and we usually work up from a base case rather than down to it. Thinking recursively, in comparison, can seem like trying to build a pyramid by first placing the capstone and then filling in stones until you reach the ground. Since the course unfortunately targets non-programmers, the Professor felt it was best to emphasize techniques which in turn emphasized the computer's role as a device which simply performs for us serial commands as we would, were we to do a problem by hand. I think it might be this sort of bias which makes adopting functional programming difficult for many programmers. I see the point, but I still think we should have taught them the recursive solution. Sometimes a bit of Magical Masonry can be fun.
(define (div-dif xs fs i j) (if (= 1 (- j i)) (/ (- (nth j fs) (nth i fs)) (- (nth j xs) (nth i xs))) (/ (- (div-dif (+ i 1) j) (div-dif i (- j 1))) (- (nth j xs) (nth i xs)))))