#! /usr/bin/env python # -*- coding: utf-8 -*- """ Python solution for the Problem IV of the CS101 Final Exam (about polynomial interpolation). This file contains some nice animations, to really understand what is happening. @date: Thu Apr 23 19:08:05 2015. @author: Lilian Besson for CS101 course at Mahindra Ecole Centrale 2015. @licence: MIT Licence (http://lbesson.mit-license.org). """ # %% Question 1) def divided_differences(X, Y, n): """ Returns a list D of size n+1, containing the values of the successive divided differences for the x points X = [x0,..,xn] and y points Y = [y0,..yn]. - Time complexity is O(2^n). - Memory complexity is O(2^n). - This function is recursive, but we could also use dynamic programming (cf. Neville's algorithm: https://en.wikipedia.org/wiki/Neville%27s_algorithm#The_algorithm). """ assert n+1 == len(X) == len(Y) if n == 0: # Base case k = 1 return [Y[0]] elif n == 1: return[Y[0], (Y[0] - Y[1])/(X[0] - X[1])] else: D1 = divided_differences(X[:-1], Y[:-1], n-1) D2 = divided_differences(X[1:], Y[1:], n-1) return D1 + [ (D1[-1] - D2[-1]) / (X[0] - X[-1]) ] # %% Question 2) def prod(iterator): """ Compute the product of the values in the iterator. - Empty product is 1. """ p = 1 for x in iterator: p *= x return p def interpolate(X, f, t): """ Compute the value of the interpolation polynomial of order $n+1$ for $f$ on the points X = [x0, .., xn], ie the value $L_{f, n}(t)$. - Time and memory complexity is O(2^n). """ n = len(X) - 1 Y = [ f(xi) for xi in X ] D = divided_differences(X, Y, n) return sum( D[i] * prod(t - X[j] for j in xrange(0, i)) for i in xrange(0, n+1) ) def L(X, f): """ L(X, f)(t) == interpolate(X, f, t) But the coefficients of the polynomials are cached so they are not computed twice. - After caching D and X, the computational time is O(n^2). - Time and memory complexity for computing D is O(2^n). """ n = len(X) - 1 Y = [ f(xi) for xi in X ] D = divided_differences(X, Y, n) return lambda t: sum( D[i] * prod(t - X[j] for j in xrange(0, i)) for i in xrange(0, n+1) ) # %% Testing all this on one example if __name__ == '__main__': print "Two examples of polynomial interpolation for sin(x)." import numpy as np import matplotlib.pyplot as plt def f(x): return np.sin(x) # n = 10 # 1000 loops, best of 3: 824 µs per loop # n = 20 # 1 loops, best of 3: 834 ms per loop # Going from n = 10 to n = 20 multiply the time by 1000 # Which is 2**10 : from 2**10 to 2**20 # So our function divided_differences is indeed in O(2^n) # n = 40 # Too big! # %timeit D = divided_differences(np.linspace(-4, 4, n), f(np.linspace(-4, 4, n)), n-1) xmin, xmax = -4, 4 n = 7 X = np.linspace(xmin, xmax, n) print "Points X = [x0, .., xn]:", X Y = f(X) print "Points Y = [y0, .., yn]:", Y D = divided_differences(X, Y, len(X)-1) print "Divided differences:", D plt.figure(1) plt.title(r"Example of polynomial interpolation for $t \mapsto \sin(t)$ with $n = {}$ points.".format(n)) # The points used for the interpolation plt.scatter(X, f(X), color='red', linewidth=6) # The function T = np.linspace(2*xmin, 2*xmax, 1000) Yf = f(T) plt.plot(T, Yf, color='red', label=r'$\sin(t)$') # The interpolation polynomial P = L(X, f) YP = [ P(t) for t in T ] plt.plot(T, YP, color='blue', label=r'$L_{f, %i}(t)$' % n) T = np.linspace(xmin, xmax, 400) error = np.sum((f(T) - np.array([ P(t) for t in T ]))**2) print "Norm 2 of the error is about {:g}.".format(error) plt.ylim(1.2*Yf.min(), 1.2*Yf.max()) plt.legend() plt.show() # %% Animation example from matplotlib import animation def makeanimation(f, xmin=0, xmax=1, namef=r'$f(x)$', frames=16): """ Make an animation for f (of name namef), from xmin to xmax.""" plt.show() fig, ax = plt.subplots() ax.set_title(r"Example of polynomial interpolation for $t \mapsto$ {} on [{}, {}]".format(namef, xmin, xmax)) ax.set_xlim(2*xmin, 2*xmax) T = np.linspace(2*xmin, 2*xmax, 1000) Tsmall = np.linspace(xmin, xmax, 400) Y0 = np.zeros_like(T) line1, = ax.plot(T, Y0, 'ro', markersize=6) line2, = ax.plot(T, Y0, linewidth=2) # Initialization function: plot the background of each frame def init(): print "init() called." # Start by plotting the function Yf = f(T) ymin, ymax = Yf.min(), Yf.max() ymin = min(1.2*ymin, 0.9*ymin) ymax = max(1.2*ymax, 0.9*ymax) ax.set_ylim(ymin, ymax) ax.plot(T, Yf, color='black', label=namef, linewidth=4) plt.vlines([xmin, xmax], [ymin, ymin], [ymax, ymax], linewidth=2, linestyle=':', color='green') # Put useless data for line1 and line2 line1.set_data([], []) # line1.set_color('red') # line1.set_linestyle('o') # line1.set_linewidth(6) line1.set_label("Interpolation points") line2.set_data([], []) # line2.set_linewidth(2) line2.set_label("Interpolation polynomial ($L_{f, n}(t)$)") ax.legend() # return line1, line2, return line2, # Animation function. This is called sequentially def animate(i): print "animate({}) called.".format(i) n = i + 1 print n # Either we take regularly spaced points: # X = np.linspace(xmin, xmax, n) # Or we take random points: X = xmin + (xmax - xmin) * np.random.random(n) print "Points X = [x0, .., xn]:", X Y = f(X) print "Points Y = [y0, .., yn]:", Y # The points used for the interpolation line1.set_color('red') line1.set_data(X, Y) line1.set_markersize(10) # The interpolation polynomial P = L(X, f) YP = [ P(t) for t in T ] line2.set_label("$L_{f, %i}(t)$" % n) line2.set_data(T, YP) ax.legend() Ysmall = f(Tsmall) YPsmall = np.array([ P(t) for t in Tsmall ]) error = np.sum(( Ysmall - YPsmall )**2) print "For {} points, estimated norm 2 error is about {:g}.".format(n, error) ax.set_title("Polynomial interpolation for $t \\mapsto$ {} on $[{}, {}]$,\n with {} points (norm 2 error is about ${:g}$).".format(namef, xmin, xmax, n, error)) if n == frames: ax.fill_between(Tsmall, Ysmall, YPsmall, color='magenta', alpha=.25) # return line1, line2, return line2, # End of animate() # Call the animator. print "Launching the animation, with {} frames.".format(frames) anim = animation.FuncAnimation( fig, animate, init_func=init, frames=frames, # FIXME: choose the best value interval=1000, # FIXME: choose the best value repeat=False, # blit=True # Only draw what have changed ) # %% Conclude by saving the animation # FIXME not working on Windows! # anim.save('Animation__Successive_polynomial_interpolation_for_ln1+x2+cosx.mp4', # fps=1, # writer='ffmpeg' # ) return anim # End of the function makeanimation # %% First test makeanimation(np.sin, xmin=-4, xmax=4, namef=r'$\sin(x)$', frames=10) # %% Second test # makeanimation(np.exp, xmin=-3, xmax=1, namef=r'$\exp(x)$', frames=5) # %% Third test # makeanimation(lambda x: np.log(1+x**2), xmin=-3, xmax=3, namef=r'$\ln(1+x^2)$', frames=15) # %% Forth test # makeanimation(lambda x: np.log(1+x**2)+np.cos(x), xmin=-5, xmax=5, namef=r'$\ln(1+x^2)+\cos(x)$', frames=20) # End