# CSE 330 Numerical Analysis Lab 

### Lab 6: Polynomial Root Finding

### Instructions
---

In this week's assignment we will try to find out roots of a polynomial.
1. Finding root via binomial method.
2. Finding root via Fixed point (Newton method).

### `pandas` Dataframe:
Before we start we will know a little about `pandas` dataframe. `pandas` is a python library. Dataframes are classes that are used to store complex data. You can initialize one as a python dictionary. Let's take a look. 


```
import pandas as pd

x = [0, 1, 2, 3, 4, 5, 6, 7]
y = [1.0, 5.2, 3.9, 8.4, 14.6, 12.3, 8.9, 10.2]

dictionary = {
    "x": x,
    "y": y,
}
df = pd.DataFrame(dictionary)
print(df)
```

       x     y
    0  0   1.0
    1  1   5.2
    2  2   3.9
    3  3   8.4
    4  4  14.6
    5  5  12.3
    6  6   8.9
    7  7  10.2


We can use it to print data in a tabular format. We can even use more `numpy` arrays with it.


```
import numpy as np

x = np.linspace(0, np.pi/2, 15)
y = np.sin(x)

dictionary = {
    "x": x,
    "sin(x)": y,
}
df = pd.DataFrame(dictionary)
print(df)
```

               x    sin(x)
    0   0.000000  0.000000
    1   0.112200  0.111964
    2   0.224399  0.222521
    3   0.336599  0.330279
    4   0.448799  0.433884
    5   0.560999  0.532032
    6   0.673198  0.623490
    7   0.785398  0.707107
    8   0.897598  0.781831
    9   1.009798  0.846724
    10  1.121997  0.900969
    11  1.234197  0.943883
    12  1.346397  0.974928
    13  1.458597  0.993712
    14  1.570796  1.000000


### Root Finding:
Let $f(x)$ be a function of $x$. if for some $x=\alpha, f(x) = f(\alpha) = 0$, we say $\alpha$ is a root of function $x$.

Let, 
$$f(x) = x^5 + 2.5x^4 - 2x^3 -6x^2 + x + 2\tag{6.1}$$

The graph of $f(x)$ looks like this.


```
from numpy.polynomial import Polynomial
import matplotlib.pyplot as plt

f = Polynomial([2.0, 1.0, -6.0, -2.0, 2.5, 1.0])
x = np.linspace(-2.5, 1.6, 100)
y = f(x)
dictionary = {
    'x': x,
    'y': y
}
plt.axhline(y=0, color='k')
plt.plot(x, y)
plt.plot(f.roots(), f(f.roots()), 'ro')
print(f.roots())
```

Roots are the points where the graph intersects with the $X$-axis.

### Bisection

One way to find out root's are to use bisection method. Here is the strategy, if $\alpha$ is a root between and interval $[a, b]$ then graph will cross the $X$-axis at $\alpha$. So, $sign( f(\alpha-h) ) = - sign( f(\alpha+h) )$, for small value of $h$. 

So, we can work our way up towards the root by taking average of $a$ and $b$, as long as the signs are different.

we will start with $a_0$ and $b_0$, such that, $f(a_0) f(b_0) < 0$.
Then we iterate as this,
\begin{align}
m_k &= \frac{a_k + b_k}{2} \\
\text{if, } f(a_k) f(m_k) < 0, \text{ then, } a_{k+1} &= a_k \text{ and } b_{k+1} = m_k\\
\text{else, } a_{k+1} &= m_k \text{ and, } b_{k+1} = b_k
\end{align}

We keep iterating until we find the root with sufficient precision. We usually use a formula like this,
$$
\frac{|m_{k+1} - m_k|}{|m_{k+1}|} \leq \epsilon \\  
$$
Where, $\epsilon$ is a very small value, like $\epsilon < 10^{-6}$

### A. Complete the code below
Complete the code below to iterate and solve for a root of $Eq. 6.1$, between the interval, $[-0.5, 1.3]$.


```
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from numpy.polynomial import Polynomial

f = Polynomial([2.0, 1.0, -6.0, -2.0, 2.5, 1.0])
a = -0.5
b = 1.3
m = (a+b)/2
e = 1e-6

roots_in_range = list(filter(lambda x: x >= a and x <= b, f.roots()))
print(f'Actual root(s): {roots_in_range} in range [{a}, {b}]')

list_a = []
list_f_a = []
list_b = []
list_f_b = []
list_m = []
list_f_m = []

root = None

# +----------------------+
# | Start of your code   +
# +----------------------+

epsilon = e + 1  # tricking the machine to start the iteration
prev_m = None  # backup variable to store previus m value
# in order to calculate epsilon later

while True:
    if prev_m is not None:
        if abs((m-prev_m)/m) >= e:
            epsilon = abs((m-prev_m)/m)  # prev_m = list_m[k-1]
        else:
            break
    prev_m = m

    list_a.append(a)  # a = list_a[k]
    list_b.append(b)  # b = list_b[k]
    list_m.append(m)  # m = list_m[k]

    f_a = f(a)
    f_b = f(b)
    f_m = f(m)

    list_f_a.append(f_a)
    list_f_b.append(f_b)
    list_f_m.append(f_m)  # f_m = listf[k]

    if f_a * f_m < 0:
        b = m
    else:
        a = m
    m = (a + b) / 2

root = prev_m  # current m value is overfitting

# +--------------------+
# | End of your code   +
# +--------------------+

print()

#  print(pd.DataFrame({'a':list_a, 'b':list_b, 'm':list_m, 'f(x)':list_f_m}))
df = {'a': list_a, 'f(a)': list_f_a,
      'b': list_b, 'f(b)': list_f_b,
      'm': list_m, 'f(m)': list_f_m}
print(pd.DataFrame(df))
print(f'\nRoot:\t\t{root}\nEpsilon:\t{epsilon}')

xs = np.linspace(-2.5, 1.6, 100)
ys = f(xs)

plt.axhline(y=0, color='k', linewidth=1.2)
plt.plot(xs, ys)
plt.plot(root, f(root), 'ro')
plt.show()
```

### Fixed point iteration

#### Fixed point
A number $\xi$ is called a **fixed point** to function $g(x)$ if $g(\xi) = \xi$. Using fixed points are a nice strategy to find roots of an equation. In this method if we are trying to find a root of $f(x) = 0$, we try to write the function in the form, $x = g(x)$. That is,
$$
f(x) = x - g(x) = 0
$$

So, if $\xi$ is a fixed point of $g(x)$ it would also be a root of $f(x)=0$, because,
$$
f(\xi) = \xi - g(\xi) = \xi - \xi = 0
$$

We can find a suitable $g(x)$ in any number of ways. Not all of them would converge, whereas, Some would converge very fast. For example, consider $Eq. 6.1$.

\begin{align}
&         &    f(x) &=x^5 + 2.5x^4 - 2x^3 -6x^2 + x + 2 \\
&\implies &x - g(x) &=x^5 + 2.5x^4 - 2x^3 -6x^2 + x + 2 \\
&\implies &    g(x) &=-x^5 - 2.5x^4 + 2x^3 + 6x^2 - 2 \tag{6.2}\\
\end{align}

again, 
$$
f(x) = x^5 + 2.5x^4 - 2x^3 -6x^2 + x + 2 = 0\\
$$
\begin{align}
&\implies &6x^2 &= x^5 + 2.5x^4 - 2x^3 + x + 2 \\
&\implies &x^2  &= \frac{1}{6}(x^5 + 2.5x^4 - 2x^3 + x + 2)\\
&\implies &x    &= \sqrt{\frac{1}{6}(x^5 + 2.5x^4 - 2x^3 + x + 2)}\\
&\implies &g(x) &= \sqrt{\frac{1}{6}(x^5 + 2.5x^4 - 2x^3 + x + 2)}\tag{6.3}\\
\end{align}

Similarly,
\begin{align}
&         &2.5x^4 &= -x^5 + 2x^3 + 6x^2 - x - 2 \\
&\implies &x^4  &= \frac{1}{2.5}(-x^5 + 2x^3 + 6x^2 - x - 2)\\
&\implies &x    &= \sqrt[\leftroot{-2}\uproot{2}4]{\frac{1}{2.5}(-x^5 + 2x^3 + 6x^2 - x - 2)}\\
&\implies &g(x) &= \sqrt[\leftroot{-2}\uproot{2}4]{\frac{1}{2.5}(-x^5 + 2x^3 + 6x^2 - x - 2)}\tag{6.4}\\
\end{align}

### B. Complete the code below
For this example we will use a couple of $g(x)$ function to find out which one converges faster.


```
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from numpy.polynomial import Polynomial

f = Polynomial([2.0, 1.0, -6.0, -2.0, 2.5, 1.0])
g1 = Polynomial([-2.0, 0.0, 6.0, 2.0, -2.5, -1.0])

def g2(x):
    p = Polynomial([2.0, 1.0, 0.0, -2.0, 2.5, 1.0])
    p_x = p(x)
#    print(p_x)
    return np.sign(p_x) * np.sqrt(np.abs(p_x)/6)

def g3(x):
    p = Polynomial([-2.0, -1.0, 6.0, 2.0, 0.0, -1.0])
    p_x = p(x)
#    print(p_x)
    return np.sign(p_x) * np.power(np.abs(p_x)/2.5, 1.0/4.0)

a1 = 0.8
g1_a = []

a2 = 0.8
g2_a = []

a3 = 0.8
g3_a = []

#+----------------------+
#| Start of your code   +
#+----------------------+

tolerance = 1e-7
for i in range(50):
    g1_a1 = g1(a1)
    g2_a2 = g2(a2)
    g3_a3 = g3(a3)
    if i > 0:
        is_repeating_value = g1_a1 - g1_a[i-1] <= tolerance
        is_repeating_value = is_repeating_value and g2_a2 - g2_a[i-1] <= tolerance
        is_repeating_value = is_repeating_value and g3_a3 - g3_a[i-1] <= tolerance
        if is_repeating_value:
            break
    a1 = g1_a1
    g1_a.append(a1)
    a2 = g2_a2
    g2_a.append(a2)
    a3 = g3_a3
    g3_a.append(a3)

print()


#+--------------------+
#| End of your code   +
#+--------------------+

xs = np.linspace(-2.5, 1.6, 100)
ys = f(xs)
dictionary = {
    'x': xs,
    'y': ys
}
plt.axhline(y=0, color='k')
plt.plot(xs, f(xs), label='f(x)')
plt.plot(xs, g1(xs), label='g1(x)')
plt.plot(xs, g2(xs), label='g2(x)')
plt.plot(xs, g3(xs), label='g3(x)')
plt.legend()
if len(g1_a) > 0:
    root = np.array([g1_a[len(g1_a)-1], g2_a[len(g2_a)-1], g3_a[len(g3_a)-1]])
    plt.plot(root, f(root), 'ro')

print(f'Roots: {f.roots()}')
print(pd.DataFrame({'g1(x)':g1_a, 'g2(x)':g2_a, 'g3(x))':g3_a,}))
```


```

```

2903072ac6f334f3e9f54e1b5209dee9a31b5ecb
ipynb
Jupyter Notebook
polynom_root.ipynb
sheikhmishar/Numerical-Analysis-Python
03a737ba38b372fb52ad773f52cd029f7da2b307
MIT
polynom_root.ipynb
sheikhmishar/Numerical-Analysis-Python
03a737ba38b372fb52ad773f52cd029f7da2b307
MIT
polynom_root.ipynb
sheikhmishar/Numerical-Analysis-Python
03a737ba38b372fb52ad773f52cd029f7da2b307
MIT
Qwen/Qwen-72B
1. YES
2. YES