Clean up the divided_difference function.
[octave.git] / divided_difference.m
1 function dd = divided_difference(f, xs)
2 ## Compute divided difference of `f` at points `xs`. The argument `xs`
3 ## is assumed to be a vector containing at least one element. If it
4 ## contains n elements, the (n-1)st divided difference will be
5 ## calculated.
6 ##
7 ## INPUTS:
8 ##
9 ## * ``f`` - The function whose divided differences we want.
10 ##
11 ## * ``xs`` - A vector containing x-coordinates. The length of `xs`
12 ## determines the order of the divided difference.
13 ##
14 ##
15 ## OUTPUTS:
16 ##
17 ## * ``dd`` - The divided difference f[xs(1), xs(2),...]
18 ##
19 order = length(xs) - 1;
20
21 if (order < 0)
22 ## Can't do anything here. Return nothing.
23 dd = NA;
24 elseif (order == 0)
25 ## Our base case.
26 dd = f(xs(1))
27 else
28 ## Order >= 1, recurse.
29
30 ## f[x0,...,x_n-1]
31 f0 = divided_difference(f, xs(1:end-1));
32 ## f[x1,...,x_n]
33 f1 = divided_difference(f, xs(2:end));
34
35 # http://mathworld.wolfram.com/DividedDifference.html
36 dd = (f0 - f1)/(xs(1) - xs(end));
37 end
38 end