## is assumed to be a vector containing at least one element. If it
## contains n elements, the (n-1)st divided difference will be
## calculated.
- order = length(xs);
-
- if (order < 1)
+ ##
+ ## INPUTS:
+ ##
+ ## * ``f`` - The function whose divided differences we want.
+ ##
+ ## * ``xs`` - A vector containing x-coordinates. The length of `xs`
+ ## determines the order of the divided difference.
+ ##
+ ##
+ ## OUTPUTS:
+ ##
+ ## * ``dd`` - The divided difference f[xs(1), xs(2),...]
+ ##
+
+ order = length(xs) - 1;
+
+ if (order < 0)
## Can't do anything here. Return nothing.
dd = NA;
- elseif (order == 1)
+ elseif (order == 0)
## Our base case.
- dd = f(xs(1))
+ dd = f(xs(1));
else
- ## Order > 1, recurse.
-
- ## f[x0,...,x_n-1]
- f0 = divided_difference(f, xs(1:end-1));
- ## f[x1,...,x_n]
- f1 = divided_difference(f, xs(2:end));
-
- # http://mathworld.wolfram.com/DividedDifference.html
- dd = (f0 - f1)/(xs(1) - xs(end))
+ ## Order >= 1.
+ cs = divided_difference_coefficients(xs);
+ dd = dot(cs, f(xs));
end
end