#! /usr/bin/env python # -*- coding: utf-8 -*- """ Python solution for the Problem IV of the CS101 Final Exam (about polynomial interpolation). For more details, and examples with plots, and animation, see the other file: CS101_Final_Exam__Problem_IV_with_animations.py @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: # We write a naive recursive function # Left divided difference: f_{n-1}[x_0, .., x_{n-1}], ..., f_{n-1}[x_0, .., x_{n-1}] D1 = divided_differences(X[0:n-1], Y[0:n-1], n-1) # Right divided difference: f_{0}[x_1], ..., f_{n-1}[x_1, .., x_n] D2 = divided_differences(X[1:n], Y[1:n], n-1) # Now we can compute the new divided difference by using its definition fn = (D1[n-1] - D2[n-1]) / (X[0] - X[n-1]) # X[-1] is the same as X[n-1] # We concatene the left difference list with the new one return D1 + [ fn ] # %% Question 2) # This has been seen in class and in labs, many times 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) # We use the maths formula, directly. return sum( D[i] * prod(t - X[j] for j in xrange(0, i)) for i in xrange(0, n+1) )