## 03 December 2008

### Interpolation in OCaml

I just wrote some code to do polynomial interpolation in O'Caml. The full darcs repository is at http://web.mit.edu/farr/www/onumerics/; here's the guts of the code. Eventually (depending on my time and level of continued interest), this may turn into a more general O'Caml numerical library. As usual, it's released under the GPL.

(*  interp.ml: Interpolation.
Copyright (C) 2008 Will M. Farr

This program is free software: you can redistribute it and/or modify
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see .
*)

let interp xs ys x =
let n = Array.length xs in
assert(Array.length ys = n);
let ytab = Array.make n ys in
for i = 1 to n - 1 do
let m = n - i in
let last_col = ytab.(i-1) and
this_col = Array.make m 0.0 in
for j = 0 to m - 1 do
let xlow = xs.(j) and
xhigh = xs.(j+i) and
ylow = last_col.(j) and
yhigh = last_col.(j+1) in
this_col.(j) <- (x-.xlow)*.yhigh/.(xhigh-.xlow) +. (x-.xhigh)*.ylow/.(xlow-.xhigh)
done;
ytab.(i) <- this_col
done;
ytab.(n-1).(0)

Here's how it works. Initially, we have two arrays of points, xs and ys, and a point at which we want to evaluate the interpolated polynomial, x. This means we have an array of 0-order interpolating polynomials (the constants ys). We can use a linear interpolation between neighboring points to construct the values of first-order interpolating polynomials at x. Then we can linearly interpolate between the neighboring first-order polynomials to construct second-order interpolating polynomials. We iterate this process until we have a single value, which is the value of the n-point interpolating polynomial at x. In the case of four points, here is a graphical representation of the procedure:

x0: y0
y01
x1: y1     y012
y12      y0123
x2: y2     y123
y23
x3: y3
The y's subscripts indicate the points involved in the interpolation. Each y is a linear interpolation of the y values above and below it to the left at the corresponding x points.