<h1>INTERPOLATIONS</h1>

<b>Group 8</b> <br>
Gardyan Priangga Akbar (2301902296)

<h2>When do we need interpolation?</h2>

Interpolation is drawing conclusions from within a set of known information. For example, if we know that 0 is the lowest number and 10 being the maximum, we can determine that the number 5 must lie in between. Interpolation has many real-life applications, such as: <br>
<ul>
    <li>When you have the cost of catering for 25 and 100 people, but you need an estimate for the cost of catering for 50 people.</li>
    <li>When deciding what laptop to buy and you know the price tag and capabilities of laptops at both the lower and higher ends, interpolation can be used to get the most optimal price and specs out of your budget.</li>
    <li>Finding the amount of employees needed to complete a task with the most optimal cost.</li>
    <li>And more...</li>
</ul>

There are a number of ways or methods to do interpolation. Two of them, for example, is Lagrange's method and Newton's Divided Difference method. In this notebook, we will be learning about these two and how we can implement them using Python. Strap youselves in, because at the end of every section there will be a playground area for you to explore and mess around! Let's go!

<h2>Lagrange Interpolation</h2>

<h3>The Theory</h3>

Lagrange's method is one of the ways for data interpolation from a set of known data points. With this method, we can interpolate the value of f(x) from any value of x from within the data set. Here is the formula:



Where: <br>
<b>n</b> = the degree of polynomial (for linear n = 1, quadratic n = 2, and so on) <br>
<b>Li(x)</b> = the weighting function

To get the weighting function, the formula is:



For some people this formula might seem quite daunting or scary even. However, this formula is just the equivalent of



<h3>Doing it in Python</h3>

First let's make a list of the data points we know.


```python
xy_values = []

#Initialize x and y values (make sure the X values are in order)
xy_values.append([0, 0])
xy_values.append([10, 227.04])
xy_values.append([15, 362.78])
xy_values.append([20, 517.35])
xy_values.append([22.5, 602.97])
xy_values.append([30, 901.67])

xy_values
```




    [[0, 0],
     [10, 227.04],
     [15, 362.78],
     [20, 517.35],
     [22.5, 602.97],
     [30, 901.67]]



Next let's decide on the order of polynomial to interpolate our data with. We will store it in a variable called <i>n</i>. For reference, to do a linear interpolation, we put our <i>n</i> value as 1. For quadratic <i>n</i> = 2, cubic <i>n</i> = 3, and so on.


```python
n = 1
```

Now let's choose a value of <i>x</i> to interpolate. Obviously, the value of <i>x</i> needs to be within our known data points, otherwise we won't be able to interpolate (that would be extrapolation).


```python
xVal = 16
```

Next we need to pick <b>two</b> points from our known data points that sandwhiches our <i>xVal</i>. We will be keeping track on the indexes. So if our <i>xVal</i> is <b>16</b>, we will be picking the x values <b>15</b> and <b>20</b> because 16 lies between them. As we see in our <i>xy_values</i> list, 15 and 20 are positioned in the indexes <b>2</b> and <b>3</b> respectively. Hence, we take a note of that in a new list.


```python
def get_first2_indexes(xy_values, xVal):
    indexes = []
    for i in range(len(xy_values)-1):
        if xy_values[i][0] < xVal and xy_values[i+1][0] > xVal:
            indexes.append(i)
            indexes.append(i+1) 
    return indexes
        
indexes = get_first2_indexes(xy_values, xVal)
indexes
```




    [2, 3]



If <i>n</i> = 1 (linear), we can go directly to finding the weighting function. However, when <i>n</i> > 1, we have to also select adjacent x values from our two chosen data points. Take note to always pick the data point closest to <i>xVal</i>.

For example when <b><i>n</i> = 3</b>:
1. Compare <b>10</b> and <b>22.5</b>
2. <b>10</b> is closer to <b>16</b> than 22.5. So we choose that.
3. <b><i>indexes</i></b> will now house [1, 2, 3]. Take note to keep track the indexes in ascending order.

For example when <b><i>n</i> = 4</b>:
1. We add one more data point from when <i>n</i> = 3.
2. Compare <b>0</b> and <b>22.5</b>
2. <b>22.5</b> is closer to <b>16</b> than 0. So we choose that.
3. <b><i>indexes</i></b> will now house [1, 2, 3, 4].


```python
def get_remaining_indexes(xy_values, indexes, xVal, n):
    for _ in range(n-1):
        #find the value nearest to xVal
        leftIndex = indexes[0]-1
        rightIndex = indexes[len(indexes)-1] + 1
        #Check if the adjacent index exists in the given xy_values data
        if (leftIndex > -1):
            if (rightIndex < len(xy_values)):
                #Check which one is closer to xVal
                if (abs(xy_values[leftIndex][0] - xVal) < abs(xy_values[rightIndex][0] - xVal)):
                    indexes.insert(0, leftIndex)
                else:
                    indexes.append(rightIndex)
            else:
                indexes.insert(0, leftIndex)
        elif (rightIndex < len(xy_values)):
            indexes.append(rightIndex)

get_remaining_indexes(xy_values, indexes, xVal, n)
            
indexes
```




    [2, 3]



Now we can go ahead and try to find the weighting functions. We will be using <b>Sympy</b> to help us keep track of variables and automatically calculate the final result. Let's start by importing Sympy library.


```python
import sympy as sp

x = sp.Symbol('x');
```

We will now proceed in determining the weighting function. Recall that the formula is




```python
def gather_weighting_functions(polynomial):
    wFunc = []      #Collection of Ln(x)
    for i in range(polynomial+1):
        subFunc = []    #Collection of individual (x - xj)/(xi-xj)
        for j in range(polynomial+1):
            #j != i
            if i != j:
                #(t - xj)/(xi-xj)
                #sub = [i, j]
                #sub[0] = xi
                #sub[1] = xj
                sub = []
                sub.append(i)
                sub.append(j)
                subFunc.append(sub)
        wFunc.append(subFunc)
        
    return wFunc

wFunc = gather_weighting_functions(n)
wFunc
```




    [[[0, 1]], [[1, 0]]]



The code above simply stores the values i and j in each of their respective iterations.

Recall the formula for lagrange's interpolation to be



We will now put <b>fn(x)</b> together with the code below (Sympy has the benefit of automatically simplifying our otherwise very long equation):


```python
def get_equation(xy_values, wFunc, indexes, x_symbol):
    total = 0
    for i in range(len(wFunc)):
        weight_function_prod = 1
        for a in range(len(wFunc[i])):
            iIndex = wFunc[i][a][0]
            index = indexes[iIndex]
            xi = xy_values[index][0]
            jIndex = wFunc[i][a][1]
            index = indexes[jIndex]
            xj = xy_values[index][0]
            sub = (x_symbol - xj)/ (xi - xj)
            weight_function_prod *= sub
        #Multiply by f(i)
        total += weight_function_prod * xy_values[indexes[i]][1]

    return sp.simplify(total)
    
equation = get_equation(xy_values, wFunc, indexes, x)
equation
```




$\displaystyle 30.914 x - 100.93$



We are not done however, because we are interested in the value of <b>y</b> when x is our <b><i>xVal</i></b>, which is <b>16</b>. To solve this, we can call Sympy's <b>evalf()</b> function on our <i>equation</i> variable. 


```python
#Solve for xVal
result = equation.evalf(subs={x : xVal})

result
```




$\displaystyle 393.694$



If we graph our findings it will look like this


```python
def graph_lagrange(xy_values, equation, xVal, result, x_symbol):
    #Graphing
    %matplotlib inline
    import matplotlib.pyplot as plt
    #split x and y
    x_values = []
    y_values = []

    for i in range(len(xy_values)):
        x_values.append(xy_values[i][0])
        y_values.append(xy_values[i][1])

    #Generate x and y
    new_x_values = []
    new_y_values = []
    for i in range(int(min(x_values) * 100), int(max(x_values) * 100), 1):
        new_x_values.append(i/100)
        new_y_values.append(equation.evalf(subs={x_symbol:i/100}))

    plt.plot(x_values, y_values, 'o', label='data')
    plt.plot(new_x_values, new_y_values, '-', label='equation')
    plt.plot([xVal], [result], '+', label="interpolated data")
    plt.legend()
    plt.xlabel("X")
    plt.ylabel("Y")

    print("y =", equation)
    plt.show()
    
graph_lagrange(xy_values, equation, xVal, result, x)
```

Notice how our interpolated <i>xVal</i> lies within our predicted equation, but not all of our known data points. There are some that are relatively quite far away from the equation line. In order to reduce this, we need to use a <b>higher</b> level polynomial (a higher value for <b>n</b>). Do take note that the highest level of polynomial you can do is equal to the number of data points you have minus 1. This is because you do not have enough data points to use a higher degree polynomial.

<h3>Try it Yourself!</h3>

Try experimenting with Lagrange's Interpolation yourself with your own data inputs. <br>
Let's start with the data points that you know:


```python
xy_values = []

#Initialize x and y values (make sure the X values are in order)
xy_values.append([0, 0])
xy_values.append([10, 227.04])
xy_values.append([15, 362.78])
xy_values.append([20, 517.35])
xy_values.append([22.5, 602.97])
xy_values.append([30, 901.67])

xy_values
```




    [[0, 0],
     [10, 227.04],
     [15, 362.78],
     [20, 517.35],
     [22.5, 602.97],
     [30, 901.67]]



Now for what value of x do you want to find?


```python
xVal = 16
xVal
```




    16



And what order of polynomial would you like to use? <br>
Note: Be sure to set your value of <b>n</b> to be the amount of data points -1. If you have 6 data points, your max value for <b>n</b> is 5.


```python
n = 3    #Order/degree of polynomial
#n = len(xy_values) - 1  #Use this to use the highest possible degree of polynomial
n
```




    3



Your inputs are now in! (Don't change anything in the code below)


```python
import sympy as sp

x = sp.Symbol('x');

indexes = get_first2_indexes(xy_values, xVal)
get_remaining_indexes(xy_values, indexes, xVal, n)
wFunc = gather_weighting_functions(n)
equation = get_equation(xy_values, wFunc, indexes, x)
result = equation.evalf(subs={x : xVal})

equation
```




$\displaystyle 0.00543466666666609 x^{3} + 0.132040000000028 x^{2} + 21.2655333333331 x - 4.25399999999809$



The code has been baked and here is the result!


```python
result
```




$\displaystyle 392.057168000002$



Now let's see how that looks like in a graph.


```python
graph_lagrange(xy_values, equation, xVal, result, x)
```

How did your graph turn out? Were you able to line up your known data points with your equation line? Maybe try with a higher value of n, or try with an entire different data set. Play around!

<h2>Newton Interpolation</h2>

<h3>The Theory</h3>

To do interpolation with Newton's method, we use Newton's Divided Difference Polynomial (NDDP) method. With this method, we wil be able to interpolate the equation of the line using the known data points. Here is the formula:



Where: <br>
<b>n</b> = the degree of polynomial <br>
<b>a<i>n</i></b> = is the divided difference

Solving for <b>a<i>n</i></b> is quite tricky, and so that part will be discussed as we learn to solve NDDP using Python.

<h3>Solving it with Python</h3>

Similar to what we did in Lagrange's, we first make a list of the data points we know.


```python
xy_values = []

#Initialize x and y values (make sure the X values are in order)
xy_values.append([0, 0])
xy_values.append([10, 227.04])
xy_values.append([15, 362.78])
xy_values.append([20, 517.35])
xy_values.append([22.5, 602.97])
xy_values.append([30, 901.67])
```

The beauty with Newton's method is that we do not need to specify what order of polynomial we want to use. That is determined by the amount of data points we have -1. So if we have 6 data points, we will be using a 5 degree polynomial. Our next step is to create the divided difference table, so let's do that.


```python
#Initialize divided difference table
def init_table():
    table = []
    
    for _ in range(len(xy_values)):
        temp = []
        for _ in range(len(xy_values) + 1):
            temp.append(-1)
        table.append(temp)
        
    #Insert x and y values to table
    for i in range(len(xy_values)):
        table[i][0] = xy_values[i][0]
        table[i][1] = xy_values[i][1]
    
    return table

table = init_table()
table
```




    [[0, 0, -1, -1, -1, -1, -1],
     [10, 227.04, -1, -1, -1, -1, -1],
     [15, 362.78, -1, -1, -1, -1, -1],
     [20, 517.35, -1, -1, -1, -1, -1],
     [22.5, 602.97, -1, -1, -1, -1, -1],
     [30, 901.67, -1, -1, -1, -1, -1]]



We first initialize an empty table by flagging all values with -1. Then we insert our known x and y values there. Our table will look something like this.



Next we need to populate the table and fill in the remaining empty cells. Here is the way to fill in the table in general:



Try taking a look at it carefully. You will see a particular pattern. Let's go ahead and fill the values in.


```python
#Do the divided difference table
def compute_table(table):
    y_bound = 1
    for col in range(2, len(table[0])):
        for row in range(y_bound, len(table)):
            try:
                delta = (table[row][col-1] - table[y_bound-1][col-1]) / (table[row][0] - table[y_bound-1][0])
            except:
                delta = 0
            #print(table[row][col-1], '-', table[y_bound-1][col-1], "divide", table[row][0], '-', table[y_bound-1][0], '=', delta)
            table[row][col] = delta
        y_bound += 1
compute_table(table)
table
```




    [[0, 0, -1, -1, -1, -1, -1],
     [10, 227.04, 22.704, -1, -1, -1, -1],
     [15, 362.78, 24.185333333333332, 0.29626666666666635, -1, -1, -1],
     [20, 517.35, 25.8675, 0.3163499999999999, 0.004016666666666713, -1, -1],
     [22.5,
      602.97,
      26.79866666666667,
      0.3275733333333335,
      0.004174222222222287,
      6.302222222222958e-05,
      -1],
     [30,
      901.67,
      30.055666666666664,
      0.36758333333333315,
      0.0047544444444444535,
      7.377777777777409e-05,
      1.434074074072601e-06]]



Next we need to get the values for our <b>a0, a1, a2,...., an</b>. Since we already did the divided difference table, we can just "steal" the values from it. This is how you can get the values from the table:




```python
#Get an values
def get_an_values(table):
    an = []
    col = 1
    for row in range(0, len(table)):
        an.append(table[row][col])
        col += 1
    return an

an = get_an_values(table)
an
```




    [0,
     22.704,
     0.29626666666666635,
     0.004016666666666713,
     6.302222222222958e-05,
     1.434074074072601e-06]



We now have all the pieces to put our NDDP puzzle together. Let's recall the formula again.



Alright let's do it in Python now.


```python
import sympy as sp
x = sp.Symbol('x')

def get_equation_newton(an, x_symb):
    func = 0
    for a in range(len(an)):
        product = an[a]
        for i in range(a):
            product *= (x_symb - xy_values[i][0])
        func += product
    
    func = sp.simplify(func)
    return func

func = get_equation_newton(an, x)

func
```




$\displaystyle x \left(1.4340740740726 \cdot 10^{-6} x^{4} - 3.3777777777671 \cdot 10^{-5} x^{3} + 0.00356481481481208 x^{2} + 0.211538888888918 x + 20.2515666666666\right)$



Great! We now have our equation line. Let's plot it and see how it looks.


```python
%matplotlib inline

import matplotlib.pyplot as plt

def graph_newton(xy_values, func, x_symb):
    #split n and y
    x_values = []
    y_values = []
    
    for i in range(len(xy_values)):
        x_values.append(xy_values[i][0])
        y_values.append(xy_values[i][1])
    
    #Generate x and y
    new_x_values = []
    new_y_values = []
    for i in range(int(min(x_values) * 100), int(max(x_values) * 100), 1):
        new_x_values.append(i/100)
        new_y_values.append(func.evalf(subs={x_symb:i/100}))
        
    plt.plot(x_values, y_values, 'o', label='data')
    plt.plot(new_x_values, new_y_values, '-', label='equation')
    plt.legend()
    plt.xlabel("X")
    plt.ylabel("Y")
    
    plt.show()
    
def graph_newton_with_interpolation(xy_values, func, x_symb, xVal, yVal):
    #split n and y
    x_values = []
    y_values = []
    
    for i in range(len(xy_values)):
        x_values.append(xy_values[i][0])
        y_values.append(xy_values[i][1])
    
    #Generate x and y
    new_x_values = []
    new_y_values = []
    for i in range(int(min(x_values) * 100), int(max(x_values) * 100), 1):
        new_x_values.append(i/100)
        new_y_values.append(func.evalf(subs={x_symb:i/100}))
        
    plt.plot(x_values, y_values, 'o', label='data')
    plt.plot(new_x_values, new_y_values, '-', label='equation')
    plt.plot([xVal], [yVal], '+', label="interpolated data")
    plt.legend()
    plt.xlabel("X")
    plt.ylabel("Y")
    
    plt.show()   

print("y =", func)
graph_newton(xy_values, func, x)
```

Well would you look at that. The graph looks quite nicely; all the data points are inside the graph line. Much better than lagrange's method using a polynomial degree of 1. And we didn't even specify the order of polynomial with Newton's.

<h3>Try it Yourself!</h3>

Try experimenting with Newton's Divided Difference method yourself with your own datasets. <br>
Let's start with the data points that you know:


```python
xy_values = []

#Initialize x and y values (make sure the X values are in order)
xy_values.append([0, 0])
xy_values.append([10, 227.04])
xy_values.append([15, 362.78])
xy_values.append([20, 517.35])
xy_values.append([22.5, 602.97])
xy_values.append([30, 901.67])

xy_values
```




    [[0, 0],
     [10, 227.04],
     [15, 362.78],
     [20, 517.35],
     [22.5, 602.97],
     [30, 901.67]]



Your inputs are now in! Let's process it.


```python
x = sp.Symbol('x')

table = init_table()
compute_table(table)
an = get_an_values(table)
func = get_equation_newton(an, x)
```

The code has been baked and here is the line equation you get:


```python
func
```




$\displaystyle x \left(1.4340740740726 \cdot 10^{-6} x^{4} - 3.3777777777671 \cdot 10^{-5} x^{3} + 0.00356481481481208 x^{2} + 0.211538888888918 x + 20.2515666666666\right)$



Let's see how it looks in a graph


```python
print("y =", func)
graph_newton(xy_values, func, x)
```

How does your graph look? Are all the data points lined up nicely? Maybe try with a higher value of <i>n</i>, or try with an entirely different data set. Play around!

Perhaps you would like to interpolate a value of y? What value for x would you like to try?


```python
xVal = 16
xVal
```




    16




```python
yVal = func.evalf(subs={x: xVal})

yVal
```




$\displaystyle 392.070578915556$



Let's see where does that lie in the graph


```python
graph_newton_with_interpolation(xy_values, func, x, xVal, yVal)
```

Does your interpolated data lie somewhere in the graph? Mess around with more interpolation and enjoy your graph :)

<h3>Conclusion</h3>

Alright, now that we've explored both Langrange's and Newton's method to do interpolations, let's do a recap!

Interpolation is all about making an educated guess from within a set of known data. Two methods we can use to do interpolation is Lagrange's and Newton's Divided Difference method. With Lagrange's we can control the degree of polynomial we want to use to produce different accuracies and therefore control the speed at which our program is going to run. Newton's method, on the other hand, would use the highest possible degree of polynomial to produce the most accurate and precise interpolation, but that means it will also use the max amount of time to compute. So which one to use? It all depends on you. Use the method that suits you best.

05e284f868b813b89b435e979697f4b6b1583279
ipynb
Jupyter Notebook
Interpolation.ipynb
GiantSweetroll/Computational-Math-Interpolation
1945649b1be4a9814deed6ab81ebd510841cd11b
Apache-2.0
Interpolation.ipynb
GiantSweetroll/Computational-Math-Interpolation
1945649b1be4a9814deed6ab81ebd510841cd11b
Apache-2.0
Interpolation.ipynb
GiantSweetroll/Computational-Math-Interpolation
1945649b1be4a9814deed6ab81ebd510841cd11b
Apache-2.0
Qwen/Qwen-72B
1. YES
2. YES