#!/usr/bin/env python
# coding: utf-8

# <h1><center>Séance 2, Interpolation, Systèmes Linéaires</h1><center>

# The general interpolation problem can be summarized as follows. We consider a set of pairs $(x_i, y_i)$, $i=0, \ldots, n$ as a well as a basis for the polynomial space, $\phi_0(x), \phi_1(x), \ldots, \phi_n(x)$. The result of the interpolation can then be written as a linear combination of the basis vectors
# \begin{align*}
# f(x) = \sum_{k=0}^n a_k\phi(x)
# \end{align*}
# 
# We want the interpolation to pass through the pairs $(x_i, y_i)$ so that $f(x_i) = y_i$. In other words, we want to determine $a_k$ such that 
# \begin{align*}
# f(x_i) = \sum_{k=0}^n a_k \phi_k(x_i) = y_i, \quad i=0, 1, \ldots, n
# \end{align*}
# which is a system of $n+1$ equations with $n+1$ unknowns. Using matrices, the problem can be written as 
# \begin{align*}
# \mathbf{A} = \left[\begin{array}{ccc}
# \phi_0(x_0)& \ldots & \phi_n(x_0)\\
# \phi_0(x_1)& \ldots & \phi_n(x_1)\\
# \vdots & & \\
# \phi_0(x_n)& \ldots & \phi_n(x_n)
# \end{array}\right] \left[\begin{array}{c}
# a_0\\
# a_1\\
# \vdots \\
# a_n
# \end{array}\right] = \left[\begin{array}{c}
# y_0\\
# y_1\\
# \vdots\\
# y_n
# \end{array}\right]
# \end{align*}
# 
# There exist several choices for the basis, including 
#  
#  - The monomial basis $\displaystyle \phi_k(x) = x^k$ 
#  - The Lagrange basis $\displaystyle \phi_k(x) = \prod_{\substack{j=0\\ j\neq k}}^n \left(\frac{x - x_j}{x_k -x_j}\right)$
#  - Newton basis $\phi_k(x) = \prod_{j=0}^{k-1} \left(x - x_j\right)$
# 
# where $k=0, \ldots, n$
# 
# In the case of the monomial basis, the matrix of the basis (which is known as the Vandermonde matrix) is given by
# 
# \begin{align*}
# \left[\begin{array}{ccccc}
# 1& x_0 & x_0^2 & \ldots & x_0^n\\
# 1& x_1 & x_1^2 & \ldots & x_1^n\\
# \vdots & & & & \\
# 1& x_n & x_n^2 & & x_n^n
# \end{array}\right]
# \end{align*}
# 
# In the case of the Lagrange basis, the matrix $\mathbf{A}$ turns into
# \begin{align*}
# \left[\begin{array}{cccc}
# \ell_0(x_0)& \ell_1(x_0)& \ldots & \ell_n(x_0)\\
# \ell_0(x_1)& \ell_1(x_1)& & \ell_n(x_1)\\
# \vdots & & & \\
# \ell_0(x_n)& \ell_1(x_n)& & \ell_n(x_n)
# \end{array}\right]
# \end{align*}
# 
# which, from the definition of the Lagrange polynomials, $\ell_k(x_k)=1$ and $\ell_n(x_i) = 0$ for all $i\neq k$ reduces to 
# \begin{align*}
# \left[\begin{array}{cccc}
# 1& 0 &\ldots & 0\\
# 0 & 1 & & \\
# \vdots & & \ddots & \\
# 0 & & & 1
# \end{array}\right]
# \end{align*}
# 
# and the corresponding system reduces to $a_k = y_k$ for all $k$.

# ### Exercice 1
# Find the interpolating polynomial using the monomial basis and Lagrange basis functions for the data $(-1,-6), (1,0), (2, 6)$ (hint use the function np.array and np/linalg.solve)

# In[ ]:


import numpy as np




# ### Exercise 2
# 
# Write down the Lagrange interpolation polynomial of degree $n$ for the given functions and nodes. Using Pyplot, plot the function, nodes and corresponding polynomial
# 
#  - $f(x) = \sin(x)$, $n=2, x_0 = 0, x_1 = \pi/2$
#  - $f(x) = \sin(x)$, $n=2$, $x_0=0, x_1 = \pi/6, x_2 = \pi/2$,
#  - $f(x) = \cos(x)$, $n=2$, $x_0 = 0, x_1 = \pi/3, x_2 = \pi/2$
#  - $f(x) = \tanh(x)$, $n=2$, $x_0 = 0, x_1 = \pi/3, x_2 = \pi/2$
# 

# In[ ]:





# ### Exercice 3
# 
# For each case below, plot the Lagrange polynomial on the interval $[x_0,x_n]$, over the interpolation points
# 
#  - $n=2$, $x_0 = -1, x_1 = -0.2, x_2 = 0$, $\ell_2(x)$
#  - $n=4$, $x_0 = 0, x_1 = 1, x_2 = 1.5, x_3 = 2.5, x_4 = 3$, $\ell_3(x)$
#  - $n = 20$, $x_i=i/n$ for $i=0,\ldots n$, $\ell_{10}(x)$

# ### Exercice 4. 
# 
# We now would like to be able to compute the Lagrange interpolating polynomial regardless of the number of points. Write a script that computes the Lagrange's interpolant polynomial using the Lagrange fundamental polynomials. 
# 
# - Start by defininig a function lagrange(k,x,z) which takes as input the index of the fundamental polynomial returned, x should be the array of interpolation points and z the point (or array of points) at which you want to evaluate the Lagrange polynomial.
# 

# In[ ]:


def lagrange(k,x,z):
    
    
    return lagrange_basis


# - Given the function lagrange(k,x,z), we then want to write a function interpolating_polynomial(x,y,z) that takes as input the interpolation points x and function values y as well as the points z at which you want to evaluate the interpolation polynomial and returns $p(z)$. 
# 

# In[ ]:


def interpolating_polynomial(x, y, z):
    
    
    return polynomial


# Test your function using the following nodes
# 
# x = np.array([-1., 0, 2, 3, 5])
# y = np.array([ 1., 3, 4, 3, 1])
# 
# x1 = np.array([-1., 0, 2, 3, 5, 6, 7])
# y1 = np.array([ 1., 3, 4, 3, 2, 2, 1])

# ### Exercise 5
# 
# For each function below, generate equispaced points on the $[0,1]$ interval. Then build the Lagrange interpolation polynomials of degree $2$, $3$ and $4$
# 
#  - $\sin(10x) + \cos(9x) - x^2 +2$ 
#  - $\log(5x +1 + \sin(10*x)) + x^3+2x^2 +1 + \cos(9x)$

# In[ ]:





# ### Exercise 6
# 
# 
# The Newton basis for polynomials up to degree $n$ is given by 
# \begin{align*}
# \pi_k(x) = \prod_{j=0}^{k-1} (x-x_j), \quad k=0,1,\ldots, n
# \end{align*}
# the polynomial $\pi_0(x) = \prod_{j=0}^{-1} (x-x_j)$ is interpreted as the constant polynomial $1$. 
# 
# From this, one can decompose a polynomial of degree $n$ as 
# 
# \begin{align*}
# p_n(x) &= a_0\pi_0(x) + a_1\pi_1(x) + \ldots +a_n\pi_n(x)\\
# & = a_0 + a_1(x-x_0) + a_2(x-x_0)(x-x_1) + \ldots + a_n(x-x_0)\ldots (x-x_{n-1})
# \end{align*}
# 
# The coefficients of this decomposition can then be obtained as
# 
# \begin{align*}
# p_n(x_i) = a_0 + a_1(x_i - x_0) + \ldots + a_n(x_i - x_0)\ldots (x_i - x_{n-1}) = y_i
# \end{align*}
# 
# Or similarly
# 
# 
# \begin{align*}
# \left[\begin{array}{ccccc}
# 1& 0 & 0 & \ldots & 0\\
# 1& (x_1 - x_0)& 0 & & 0\\
# 1& (x_1 - x_0) & (x_2 - x_0) (x_2 - x_1)& & 0\\
# \vdots & \vdots & \vdots & & \vdots\\
# 1& (x_n-x_0) & (x_n - x_0)(x_n-x_1) & \ldots & \prod_{i=0}^{n-1} (x_n - x_i)
# \end{array}\right]\left[\begin{array}{c}
# a_0\\
# a_1\\
# \vdots \\
# a_n
# \end{array}\right] = \left[\begin{array}{c}
# y_0\\
# y_1\\
# \vdots \\
# y_n
# \end{array}\right]
# \end{align*}
# 

# Write a script that computes the interpolating polynomial recursively using the divided differences by (a) solving the linear system and (b) computing the coefficients recursively. Then use your script to find the interpolating polynomial using the Newton basis for the data $(-1,-6), (1,0)$ and $(2,6)$. Plot this polynomial on top of the samples.

# In[ ]:


def NewtonPolynomial_linearSystem(x, y):
    
    
    
    
    return polynomial_coeffs
    
    
    
def NewtonPolynomial_Recursive(x, y):
    
    
    
    
    return polynomial_coeffs
    
    


# ### Exercise 7 Runge phenomenon
# 
# We consider the function $f(x) = \frac{1}{1+25x^2}$. Generate 15 equispaced points on the [-1,1] interval and plot the resulting Lagrange interpolation polynomial on top of the function and the sample. What do you notice ? How can we explain this phenomenon?

# In[ ]:





# Repeat the previous step with Chebyshev interpolation points. For this, write a function Chebyshev(a, b) that, given an interval $[a, b]$ returns the Chebyshev interpolation points. 

# In[ ]:


def Chebyshev(a, b):
    
    
    
    
    return ChebInterp


# ### Exercise 8
# 
# On souhaite predire la quantite de cabillaud qui sera pechee au niveau Europeen en 2026. Pour ce faire on dispose du tableau suivant (FAO) qui reprend les quantites pechees en tonne
#     
# \begin{align*}
# \begin{array}{|l|lllll|}\hline
#   Année &1950&1960&1970&1980&1990 & 2000 & 2010 & 2020\\ \hline
#     \text{Quantité déchargées ($10^6$ T)} &1.32&1.76&2.21& 1.53&1.1 & 0.8 & 0.85 & 1  \\
#     \hline
#   \end{array}
#   \end{align*}
#   
# En utilisant les années 1950 a 2010, on commencera par predire, a l'aide d'un modèle d'interpolation polynomial, le déchargement de cabillaud en 2020. On calculera ensuite l'erreur commise. On predira a l'aide du même modèle, la production en 2026. Finalement on comparera graphiquement l'evolution du dechargement exact (voir ci dessous) a l'evolution predite par le modele sur les annees suivantes:
# 
# \begin{align*}
# \begin{array}{|l|lllll|}\hline
#   Année &1950&1955&1960&1965&1970 & 1975 & 1980 & 1985& 1990&1995&2000&2005&2010 & 2015 & 2020 \\ \hline
#     \text{Q($10^6$ T)} &1.32& 1.72 & 1.76& 1.91 & 2.21& 1.68& 1.53& 1.4 &1.1 & 1.34 & 0.8 & 0.73& 0.85 & 1.32 & 1  \\
#     \hline
#   \end{array}
#   \end{align*}

# In[ ]:





# ### Exercice 9
# 
# Implement the Neville-Aitken scheme. The function should return a list of all the polynomials. You can for example use a 2D array whose first column contains the $y_i$, then second column contains the degree 1 interpolating polynomials, ...

# In[ ]:


def Neville_Aitken(x, y)
   
   
   return list_coeff