Add the divided_difference() and newton_polynomial() functions to interpolation.py.

diff --git a/mjo/interpolation.py b/mjo/interpolation.py
index e3e69338348ccdedb9fc581569f71422206cf787..747c680da2f73c777c2f86cf4f20e4cb4300b3a1 100644 (file)
--- a/mjo/interpolation.py
+++ b/mjo/interpolation.py
@@ -71,3 +71,93 @@ def lagrange_polynomial(f, x, xs):
     ls = [ lagrange_coefficient(k, x, xs) for k in range(0, len(xs)) ]
     sigma = sum([ ys[k] * ls[k] for k in range(0, len(xs)) ])
     return sigma
+
+
+def divided_difference(f, xs):
+    """
+    Return the Newton divided difference of `f` at the points
+    `xs`. Reference:
+
+      http://en.wikipedia.org/wiki/Divided_differences
+
+    INPUT:
+
+      - ``f`` -- The function whose divided difference we seek.
+
+      - ``xs`` -- The list of points at which to compute `f`.
+
+    OUTPUT:
+
+    The divided difference of `f` at ``xs``.
+
+    TESTS::
+
+        sage: divided_difference(sin, [0])
+        0
+        sage: divided_difference(sin, [0, pi])
+        0
+        sage: divided_difference(sin, [0, pi, 2*pi])
+        0
+
+    We try something entirely symbolic::
+
+        sage: f = function('f', x)
+        sage: divided_difference(f, [x])
+        f(x)
+        sage: x1,x2 = var('x1,x2')
+        sage: divided_difference(f, [x1,x2])
+        (f(x1) - f(x2))/(x1 - x2)
+
+    """
+    if (len(xs) == 1):
+        # Avoid that goddamned DeprecationWarning when we have a named
+        # argument but don't know what it is.
+        if len(f.variables()) == 0:
+            return f(xs[0])
+        else:
+            v = f.variables()[0]
+            return f({ v: xs[0] })
+
+    # Use the recursive definition.
+    numerator  = divided_difference(f, xs[1:])
+    numerator -= divided_difference(f, xs[:-1])
+    return numerator / (xs[-1] - xs[0])
+
+
+def newton_polynomial(f, x, xs):
+    """
+    Return the Newton form of the interpolating polynomial of `f` at
+    the points `xs` in the variable `x`.
+
+    INPUT:
+
+      - ``f`` -- The function to interpolate.
+
+      - ``x`` -- The independent variable to use for the interpolating
+        polynomial.
+
+      - ``xs`` -- The list of points at which to interpolate `f`.
+
+    OUTPUT:
+
+    A symbolic function.
+
+    TESTS:
+
+        sage: xs = [ -pi/2, -pi/6, 0, pi/6, pi/2 ]
+        sage: L = lagrange_polynomial(sin, x, xs)
+        sage: N = newton_polynomial(sin, x, xs)
+        sage: bool(N == L)
+        True
+
+    """
+    degree = len(xs) - 1
+
+    N = SR(0)
+
+    for k in range(0, degree+1):
+        term  = divided_difference(f, xs[:k+1])
+        term *= product([ x - xk for xk in xs[:k]])
+        N += term
+
+    return N