From: Michael Orlitzky <michael@orlitzky.com>
Date: Fri, 26 Oct 2012 01:07:12 +0000 (-0400)
Subject: Add the divided_difference() and newton_polynomial() functions to interpolation.py.
X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=0f5a866077886101e88d89eb24d8902a4a2e508a;p=sage.d.git

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

diff --git a/mjo/interpolation.py b/mjo/interpolation.py
index e3e6933..747c680 100644
--- 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