Make sure you remove `raise NotImplementedError()` and fill in any place that says `# YOUR CODE HERE`, as well as your `NAME`, `ID`, and `LAB_SECTION` below:


```python
NAME = "FAHMID BIN KIBRIA"
ID = "19201063"
LAB_SECTION = "02"
```

---

# Part 1: Hermite Interpolation
---
Hermite Interpolation is an example of a variant of the interpolation problem, where the interpolant matches one or more **derivatives of $f$**  at each of the nodes, in addition to the function values.

## Importing the necessary libraries


```python
import numpy as np
import matplotlib.pyplot as plt
from itertools import combinations
from numpy.polynomial import Polynomial
```

## Creating the components for Hermite interpolation

For the case of Hermite Interpolation, we look for a polynomial that matches both $f'(x_i)$ and $f(x_i)$ at the nodes $x_i = x_0,\dots,x_n$. Say you have $n+1$ data points, $(x_0, y_0), (x_1, y_1), x_2, y_2), \dots, (x_n, y_n)$ and you happen to know the first-order derivative at all of these points, namely, $(x_0, y_0 ^\prime ), (x_1, y_1 ^\prime ), x_2, y_2 ^\prime ), \dots ,(x_n, y_n ^\prime )$. According to hermite interpolation, since there are $2n + 2$ conditions; $n+1$ for $f(x_i)$ plus $n+1$ for $f'(x_i)$; you can fit a polynomial of order $2n+1$. 

General form of a $2n+1$ degree Hermite polynomial:

$$p_{2n+1} = \sum_{k=0}^{n} \left(f(x_k)h_k(x) + f'(x_k)\hat{h}_k(x)\right), \tag{1}$$

where $h_k$ and $\hat{h}_k$ are defined using Lagrange basis functions by the following equations:

$$h_k(x) = (1-2(x-x_k)l^\prime_k(x_k))l^2_k(x_k), \tag{2}$$

and

$$\hat{h}_k(x) = (x-x_k)l^2_k(x_k), \tag{3}$$

where the Lagrange basis function being:

$$l_k(x) = \prod_{j=0, j\neq k}^{n} \frac{x-x_j}{x_k-x_j}. \tag{4}$$

**Note** that, we can rewrite Equation $(2)$ in this way,

\begin{align}
h_k(x) &= \left(1-2(x-x_k)l^\prime_k(x_k) \right)l^2_k(x_k) \\
&= \left(1 - 2xl^\prime_k(x_k) + 2x_kl^\prime_k(x_k) \right)l^2_k(x_k) \\
&= \left(1 + 2x_kl^\prime_k(x_k) - 2l'_k(x_k)x \right) l^2_k(x_k) \tag{5}
\end{align}
Replacing $l^\prime_k(x_k)$ with $m$, we get:
$$h_k(x) = (1 - 2xm + 2x_km)l^2_k(x_k). \tag{6}$$

# Tasks:

* The functions: `l(k, x)`, `h(k, x)` and `h_hat(k, x)` calculate the corresponding $l_k$, $h_k$, and $\hat{h}_k$, respectively.

* Function `l(k, x)` has already been defined for you. Your task is to complete the `h(k, x)`, `h_hat(k, x)`, and `hermit(x, y, y_prime)` functions.

* Later we will draw some plots to check if the code is working.

---

### Part 1: Calculate $l_k$
This function uses the following equation to calculate $l_k(x)$ and returns a polynomial:

$$l_k(x) = \prod_{j=0, j\neq k}^{n} \frac{x-x_j}{x_k-x_j}.$$


```python
# Already written for you.

def l(k, x):
    n = len(x)
    assert (k < len(x))
    
    x_k = x[k]
    x_copy = np.delete(x, k)
    
    denominator = np.prod(x_copy - x_k)
    
    coeff = []
    
    for i in range(n):
        coeff.append(sum([np.prod(x) for x in combinations(x_copy, i)]) * (-1)**(i) / denominator)
    
    coeff.reverse()
    
    return Polynomial(coeff)
```

### Part 2: Calculate $h_k$
This function calculates $h_k(x)$ using the following equation:
$$h_k(x) = \left(1 + 2x_kl^\prime_k(x_k) - 2l'_k(x_k)x \right) l^2_k(x_k).$$

This equation is basically a multiplication of two polynomials.

First polynomial: $1 + 2x_kl^\prime_k(x_k) - 2l'_k(x_k)x$.

Second polynomial: $l^2_k(x_k)$.

The `coeff` variable should contain a python list of coefficient values for the **first** polynomial of the equation. These coefficient values are used to create a polynomial `p`.


```python
def h(k, x):
    # initialize with None. Replace with appropriate values/function calls
    # initialize with None. Replace with appropriate values/function calls
    l_k = None
    l_k_sqr = None
    l_k_prime = None
    coeff = None
    p = None

    # --------------------------------------------
    # YOUR CODE HERE
    l_k = l(k,x)
    l_k_sqr = l_k**2
    l_k_prime = l_k.deriv()
    coeff = [ 1 + 2*x[k]*l_k_prime(x[k]) , -2*l_k_prime(x[k])]
    p = Polynomial(coeff)
    # --------------------------------------------
    
    return p * l_k_sqr
```


```python
# Test case for the h(k, x) function
x = [3, 5, 7, 9]
k = 2
h_test = h(k, [3, 5, 7, 9])
h_result = Polynomial([-2.5, 0.5]) * (l(k, x) ** 2)
```

### Part 3: Calculate $\hat{h}_k$
This function calculates $\hat{h}_k(x)$ using the following equation:

$$\hat{h}_k(x) = (x-x_k)l^2_k(x_k).$$

This equation is also a multiplication of two polynomials.

First polynomial: $x-x_k$.

Second polynomial:  $l^2_k(x_k)$.

The `coeff` variable should contain a python list of coefficient values for the **first** polynomial of the equation. These coefficient values are used to create a polynomial `p`.


```python
def h_hat(k, x):
    # Initialize with none
    l_k = None
    l_k_sqr = None
    coeff = None
    p = None
    
    # --------------------------------------------
    # YOUR CODE HERE
    l_k = l(k,x)
    l_k_sqr = l_k**2
    coeff = [-x[k] , 1]
    p = Polynomial(coeff)

    # --------------------------------------------
    
    return p * l_k_sqr
```


```python
# Test case for the h(k, x) function
x = [3, 5, 7, 9]
k = 2
h_test = h_hat(k, [3, 5, 7, 9])
h_result = Polynomial([-7, 1]) * (l(k, x) ** 2)
```

### Part 4: The Hermite Polynomial
This function uses the following equation:

$$p_{2n+1} = \sum_{k=0}^{n} \left(f(x_k)h_k(x) + f'(x_k)\hat{h}_k(x)\right).$$

The polynomial denoted by the equation is calculated by the variable `f`.


```python
def hermit(x, y, y_prime):
    assert len(x) == len(y)
    assert len(y) == len(y_prime)
    
    f = 0
    # --------------------------------------------
    # YOUR CODE HERE
    for i in range(len(y)):
      f += y[i]*h(i,x) + y_prime[i]*h_hat(i,x)

    print(f)
    # --------------------------------------------
    return f
```

## Testing our methods by plotting graphs.

**Note:** 

* For each of the 5 plots, there will be 2 curves plotted: one being the original function, and the other being the interpolated curve. 

* The original functions are displayed in orange color, while the hermite interpolated curves are in blue.

* `x`, `y`, and `y_prime` contain $x_i$, $f(x_i)$, and $f'(x_i)$ of the given nodes of the original function $f$.

Upon calling the `hermit()` function, it returns a polynomial `f`. For example, for plot 1, it is called `f3`.

In general, a polynomial may look like the following: $f = 1 + 2x + 3x^2$. Next, we pass in a number of $x$ values to the polynomial by calling the `.linspace()` function on the polynomial object using `f.linspace()`. This function outputs a tuple, which is stored in a variable called `data`. First element of `data` contains a 1D numpy array of $x_i$ values generated by `linspace()`, and the second element of `data` contains a 1D numpy array of the corresponding $y_i$ values outputted by our example polynomial:
$f = 1 + 2x + 3x^2$. 

Using `test_x`, we generate a range of $x_i$ values to plot the original function, and `test_y` contains the corresponding $y_i$ values of the original function. For the first plot, our original function is the *sine curve*.

For all the plots:

`plt.plot(test_x, test_y)` plots the original function.

`plt.plot(data[0], data[1])` plots the interpolated polynomial.


```python
pi      = np.pi
x       = np.array([0.0, pi/2.0,  pi, 3.0*pi/2.0])
y       = np.array([0.0,    1.0, 0.0,       -1.0])
y_prime = np.array([1.0,    0.0, 1.0,        0.0])
```

**Plot 1:** trying to interpolate a sine curve (`np.sin()`) using first 2 nodes in `x` and `y`, and their corresponding derivative in `y_prime`.


```python
n      = 1
f3     = hermit(x[:(n+1)], y[:(n+1)], y_prime[:(n+1)])
data   = f3.linspace(n=50, domain=[-3, 3])
test_x = np.linspace(-3, 3, 50, endpoint=True)
test_y = np.sin(test_x)

plt.plot(data[0], data[1])
plt.plot(test_x, test_y)
plt.show()
```

**Plot 2:** trying to interpolate a sine curve (`np.sin()`) using first 3 nodes in `x` and `y` and their corresponding derivative in `y_prime`.


```python
n      = 2
f5     = hermit(x[:(n+1)], y[:(n+1)], y_prime[:(n+1)])
data   = f5.linspace(n=50, domain=[-0.7, 3])
test_x = np.linspace(-2*pi, 2*pi, 50, endpoint=True)
test_y = np.sin(test_x)

plt.plot(test_x, test_y) # 25-
plt.plot(data[0], data[1]) # 10-33
plt.show()

```

**Plot 3:** trying to interpolate a sine curve (`np.sin()`) using first 4 nodes in `x` and `y` and their corresponding derivative in `y_prime`.


```python
n      = 3
f7     = hermit(x[:(n+1)], y[:(n+1)], y_prime[:(n+1)])
data   = f7.linspace(n=50, domain=[-0.3, 3])
test_x = np.linspace(-2*pi, 2*pi, 50, endpoint=True)
test_y = np.sin(test_x)

plt.plot(data[0], data[1])
plt.plot(test_x, test_y)
plt.show()

```

**Plot 4:** trying to interpolate an exponential curve (`np.exp()`) using all nodes in `x` and `y` and their corresponding derivatives in `y_prime`.


```python
#defining new set of given node information: x, y and y'
x       = np.array([0.0, 1.0,          2.0       ])
y       = np.array([1.0, 2.71828183,  54.59815003])
y_prime = np.array([0.0, 5.43656366, 218.39260013])


f7      = hermit( x, y, y_prime)
data    = f7.linspace(n=50, domain=[-0.5, 2.2])
test_x  = np.linspace(-0.5, 2.2, 50, endpoint=True)
test_y  = np.exp(test_x**2)

plt.plot(data[0], data[1])
plt.plot(test_x, test_y)
plt.show()


```

**Plot 5:** trying to interpolate $y = (x-3)^2 + 1$ using all nodes in `x` and `y` and their corresponding derivatives in `y_prime`.

For this plot you might be able to see only one curve due to the two curves overlapping. This means that our polynomial is accurately interpolating the original function.



```python
#defining new set of given node information: x, y and y'
x       = np.array([1.0, 3.0, 5.0])
y       = np.array([5.0, 1.0, 5.0])
y_prime = np.array([-4.0, 0.0, 4.0])

f7      = hermit( x, y, y_prime)
data    = f7.linspace(n=50, domain=[-10, 10])
test_x  = np.linspace(-10, 10, 50, endpoint=True)
test_y  = (test_x-3)**2 + 1

plt.plot(data[0], data[1])
plt.plot(test_x, test_y)
plt.show()
```

## Part 2: Polynomial Interpolation Using Newton's Divided Difference Form
---


### Newton's Divided Difference Form

Newton form of a $n$ degree polynomial:

$$p_n(x) = \sum_{k=0}^{n} a_kn_k(x),$$
where the basis is:
$$n_k(x) = \prod_{j=0}^{k-1}(x-x_j),$$
$$ n_0(x)=1,$$

and the coefficients are: $$a_k = f[x_0, x_1, ..., x_k],$$

where the notation $f[x_0, x_1,\dots,x_k]$ denotes the divided difference.

By expanding the Newton form, we get:

$$p(x) = f [x_0] + (x-x_0) f[x_0,x_1] + (x-x_0) (x-x_1) f[x_0,x_1,x_2] + \dots + (x-x_0) (x-x_1) \dots (x-x_{k-1}) f[x_0, x_1, \dots, x_k]$$

### Tasks: 
1. Complete the `calc_div_diff(x,y)` function which takes input `x` and `y`, and calculates all the divided differences. You may use the lambda function `difference()` inside the `calc_div_diff(x,y)` function to calculate the divided differences.

2. Complete the `__call__()` function which takes an input `x`, and calculates `y` using all the difference coefficients. `x` can be a single value or a numpy. In this case, it is a numpy array.

`res` variable must contain all results (corresponding y for x).


```python
class Newtons_Divided_Differences:
  
    def __init__(self, differences):
        self.differences = differences

    def __call__(self, x):
        '''
        this function is for calculating y from given x using all the difference coefficients
        x can be a single value or a numpy
        the formula being used:
        f(x) = f [x0] + (x-x0) f[x0,x1] + (x-x0) (x-x1) f[x0,x1,x2] + . . . + (x-x0) (x-x1) . . . (x-xk-1) f[x0, x1, . . ., xk]

        work on this after implementing 'calc_div_diff'. Then you should have
        f[x0], f[x0,x1]. . . . . ., f[x0, x1, . . ., xk] stored in self.differences

        'res' variable must return all the results (corresponding y for x)
        '''

        res = np.zeros(len(x)) #Initialization to avoid runtime error. You can change this line if you wish

        #----------------------------------------------
        # YOUR CODE HERE
        for k in range(0, len(x)):
          res[k]=self.differences[0]
          for i in range(1, len(self.differences)):
            product = 1
            for j in range(i):  
              product = product * (x[k] - data_x[j])
            res[k] = res[k] + product * self.differences[i]
        
        
        #----------------------------------------------
        
        return res

```


```python
# basic rule for calculating the difference, implanted in the lambda function. 
# You may use it if you wish
difference = lambda y2, y1, x2, x1: (y2-y1)/(x2-x1)

def calc_div_diff(x,y):
    assert(len(x)==len(y))
    #write this function to calculate all the divided differences in the list 'b'
    b = []  #initializing
    #----------------------------------------------
    # YOUR CODE HERE
    table=np.array([[0.000000000]*len(y)]*len(y))
    for i in range(0,len(y)-1):
        table[i, 0]=y[i]
    
    for i in range(1, len(y)):  
        for j in range(len(y) - i):  
            table[j, i] = difference(table[j, i - 1], table[j + 1, i - 1], x[j], x[i + j])
    b = [0.0000000]*len(table)
    for i in range(0,len(table)):
        b[i]=table[0, i]
    #----------------------------------------------
    return b
```


```python
# Test case for the calc_div_diff(x,y) function.

data_x = [-3.,-2.,-1.,0.,1.,3.,4.]
data_y = [-60.,-80.,6.,1.,45.,30.,16.]

test = calc_div_diff(data_x, data_y)

assert len(test) == len(data_x)

```

### Plotting the polynomial
* `data_x` and `data_y` are the coordinates of the given nodes.

* `differences` is a list which contains the divided differences as each of its elements: $f[x_0], f[x_0,x_1], f[x_0,x_1,x_2], \dots$

* `obj` is an object of type `Newtons_Divided_Differences`. Creating the object runs the constructor of the class where the `difference` are stored in `self.differences`.

* `X` contains $x_i$ values through which we want to plot our polynomial.

* Calling the object using `obj(X)` executes the `__call__()` function of the class, which returns a numpy array containing the corresponding $y_i$ values, and storing them in variable `F`.

* Using `plt.plot(X,F)`, we plot the $(x_i, y_i)$ pairs of the polynomial.


```python
import numpy as np
import matplotlib.pyplot as plt

data_x = np.array([-3.,-2.,-1.,0.,1.,3.,4.])
data_y = np.array([-60.,-80.,6.,1.,45.,30.,16.])
differences = calc_div_diff(list(data_x), list(data_y))
p = Newtons_Divided_Differences(list(differences))
test_x = np.linspace(-3, 4, 50, endpoint=True)
test_y = p(test_x)

#generating 50 points from -3 to 4 in order to create a smooth line
plt.plot(test_x, test_y)
plt.plot(data_x, data_y, 'ro')
plt.show()

```


```python

```

6e9a1ac8a275590a440a3715de9bbef4cf2f592b
ipynb
Jupyter Notebook
Hermite_and_Newtons_Divided_Difference_Interpolation_Final.ipynb
Fahmid005/CSE330-Lab-
f3537ab2ef853caceb66a42645cbc742b9a8fa68
MIT
Hermite_and_Newtons_Divided_Difference_Interpolation_Final.ipynb
Fahmid005/CSE330-Lab-
f3537ab2ef853caceb66a42645cbc742b9a8fa68
MIT
Hermite_and_Newtons_Divided_Difference_Interpolation_Final.ipynb
Fahmid005/CSE330-Lab-
f3537ab2ef853caceb66a42645cbc742b9a8fa68
MIT
Qwen/Qwen-72B
1. YES
2. YES