# 13 Root Finding

An important tool in the computational tool box is to find roots of equations for which no closed form solutions exist:

We want to find the roots $x_0$ of

$$
f(x_0) = 0
$$

## Problem: Projectile range 
The equations of motion for the projectile with linear air resistance (see *12 ODE applications*) can be solved exactly.

As a reminder: the linear drag force is
$$
\mathbf{F}_1 = -b_1 \mathbf{v}\\
b := \frac{b_1}{m}
$$

Equations of motion with force due to gravity $\mathbf{g} = -g \hat{\mathbf{e}}_y$

\begin{align}
\frac{d\mathbf{r}}{dt} &= \mathbf{v}\\
\frac{d\mathbf{v}}{dt} &= - g \hat{\mathbf{e}}_y -b \mathbf{v} 
\end{align}

### Analytical solution of the equations of motion
(Following Wang Ch 3.3.2)

Solve $x$ component of the velocity 

$$
\frac{dv_x}{dt} = -b v_x
$$

by integration:

$$
v_x(t) = v_{0x} \exp(-bt)
$$

The drag force reduces the forward velocity to 0.

Integrate again to get the $x(t)$ component

$$
x(t) = x_0 + \left[1 - \frac{v_{0x}}{b} \exp(-bt)\right]
$$

Integrating the $y$ component of the velocity

$$
\frac{dv_y}{dt} = -g - b v_y
$$

gives

$$
v_y = \left(v_{0y} + \frac{g}{b}\right) \exp(-bt) - \frac{g}{b}
$$

and integrating again

$$
y(t) = y_0 + \frac{v_{0y} + \frac{g}{b}}{b} \left[1 - \exp(-bt)\right] - \frac{g}{b} t
$$

(Note: This shows immediately that the *terminal velocity* is

$$
\lim_{t\rightarrow\infty} v_y(t) = - \frac{g}{b},
$$

i.e., the force of gravity is balanced by the drag force.)

#### Analytical trajectory 

To obtain the **trajectory $y(x)$** eliminate time (and for convenience, using the origin as the initial starting point, $x_0 = 0$ and $y_0 = 0$. Solve $x(t)$ for $t$

$$
t = -\frac{1}{b} \ln \left(1 - \frac{bx}{v_{0x}}\right)
$$

and insert into $y(t)$:

$$
y(x) = \frac{x}{v_{0x}} \left( v_{0y} + \frac{g}{b} \right) + \frac{g}{b^2} \ln \left(1 - \frac{bx}{v_{0x}}\right)
$$

#### Plot 

Plot the analytical solution for $\theta = 30^\circ$ and $v_0 = 100$~m/s.


```python
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use('ggplot')
```


```python
def y_lindrag(x, v0, b1=0.2, g=9.81, m=0.5):
    b = b1/m
    v0x, v0y = v0
    return x/v0x * (v0y + g/b) + g/(b*b) * np.log(1 - b*x/v0x)

def initial_v(v, theta):
    x = np.deg2rad(theta)
    return v * np.array([np.cos(x), np.sin(x)])
```


```python
X = np.concatenate([np.linspace(0, 42, 100), np.linspace(42, 45, 1000)])
```


```python
Y = y_lindrag(X, initial_v(100, 30), b1=1)
```

    /Users/oliver/anaconda3/lib/python3.5/site-packages/ipykernel/__main__.py:4: RuntimeWarning: invalid value encountered in log



```python
plt.plot(X, Y)
```

Compare to the numerical solution (from 12 ODE Applications):


```python
import ode

def simulate(v0, h=0.01, b1=0.2, g=9.81, m=0.5):
    def f(t, y):
        # y = [x, y, vx, vy]
        return np.array([y[2], y[3], -b1/m * y[2], -g - b1/m * y[3]])

    vx, vy = v0
    t = 0
    positions = []
    y = np.array([0, 0, vx, vy], dtype=np.float64)
    
    while y[1] >= 0:
        positions.append([t, y[0], y[1]])  # record t, x and y
        y[:] = ode.rk4(y, f, t, h)
        t += h
        
    return np.array(positions)
```


```python
r = simulate(initial_v(100, 30), h=0.01, b1=1)
```


```python
plt.plot(X, Y, lw=2, label="analytical")
plt.plot(r[:, 1], r[:, 2], '--', label="RK4")
plt.legend(loc="best")
plt.xlabel("$x$ (m)"); plt.ylabel("$y$ (m)")
```

### Predict the range $R$
How far does the ball or projectile fly, i.e., that value $x=R$ where $y(R) = 0$:

$$
\frac{R}{v_{0x}} \left( v_{0y} + \frac{g}{b} \right) + \frac{g}{b^2} \ln \left(1 - \frac{bR}{v_{0x}}\right) = 0
$$

This *transcendental equation* can not be solved in terms of elementary functions.

Use a **root finding** algorithm.

## Root-finding with the Bisection algorithm
**Bisection** is the simplest (but very robust) root finding algorithm that uses trial-and-error:

* bracket the root
* refine the brackets
* see [13_Root-finding-algorithms (PDF)](13_Root-finding-algorithms.pdf)

More specifically
1. determine a bracket that contains the root: $x_{-} < x_0 < x_{+}$
2. cut bracket in half: $x' = \frac{1}{2}(x_{+} + x_{-})$
3. determine in which half the root lies: $f(x_{-}) f(x') < 0$ then the root lies in the left half, otherwise the right half.
4. repeat until $x_{+} - x_{-} < \epsilon$.

### Implementation of Bisection

- Test that the initial bracket contains a root; if not, return `None` (and possibly print a warning).
- If either of the bracket points is a root then return the bracket point.
- Allow `Nmax` iterations or until the convergence criterion `eps` is reached.
- Bonus: print a message if no root was found after `Nmax` iterations, but print the best guess and the error (but return `None`).



```python
def bisection(f, xminus, xplus, Nmax=100, eps=1e-14):
    fa, fb = f(xminus), f(xplus)
    if (fa*fb) > 0:
        print("bisect: Initial bracket [{0}, {1}] "
              "does not contain a single root".format(xminus, xplus))
        return None
    if fa == 0:
        return xminus
    if fb == 0:
        return xplus
    
    for iteration in range(Nmax):
        x = (xplus + xminus)/2
        fx = f(x)
        if f(xplus) * fx > 0:
            # root is not between xplus and x
            xplus = x
        else:
            xminus = x
        if np.abs(fx) < eps:
            break
    else:
        print("bisect: no root found after {0} iterations (eps={1}); "
              "best guess is {2} with error {3}".format(Nmax, eps, x, fx))
        x = None
    return x
```

### Finding the range with the bisection algorithm

Define the trial function:


```python
v = initial_v(100, 30)

def f(x):
    return y_lindrag(x, v, b1=b1)
```

The initial bracket is a little bit difficult for this function: choose the right bracket near the point where the argument of the logarithm becomes 0:

$$
x_{+} = \frac{v_{0x}}{b} - \epsilon
$$


```python
b1 = 1.
m = 0.5
b = b1/m
bisection(f, 0.1, v[0]/b - 1e-12, eps=1e-6)
```




    43.300674233470772



### Find the range as a function of the initial angle 


```python
b1 = 1.
m = 0.5
b = b1/m
v0 = 100
u = []
for theta in np.arange(1, 90):
    v = initial_v(v0, theta)
    def f(x):
        return y_lindrag(x, v, b1=b1)
    R = bisection(f, 0.1, v[0]/b - 1e-16, eps=1e-5)
    if R is not None:
        u.append((theta, R))
u = np.array(u)
```

    /Users/oliver/anaconda3/lib/python3.5/site-packages/ipykernel/__main__.py:4: RuntimeWarning: divide by zero encountered in log



```python
plt.plot(u[:, 0], u[:, 1])
```

Write a function `find_range()` to calculate the range for a given initial velocity $v_0$ and plot $R(\theta)$ for $10\,\text{m/s} ≤ v_0 ≤ 100\,\text{m/s}$.


```python
def find_range(v0, b1=1, m=0.5):
    b = b1/m
    u = []
    for theta in np.arange(1, 90):
        v = initial_v(v0, theta)
        def f(x):
            return y_lindrag(x, v, b1=b1)
        R = bisection(f, 0.1, v[0]/b - 1e-16, eps=1e-5)
        if R is not None:
            u.append((theta, R))
    return np.array(u)
```


```python
for v0 in (10, 25, 50, 75, 100):
    u = find_range(v0)
    plt.plot(u[:, 0], u[:, 1], label="{} m/s".format(v0))
plt.legend(loc="best")
plt.xlabel(r"$\theta$ (degrees)")
plt.ylabel(r"range $R$ (m)")
```

## Newton-Raphson algorithm
(see derivation in class and in the PDF or [Newton's Method](http://mathworld.wolfram.com/NewtonsMethod.html) on MathWorld)


### Activity: Implement Newton-Raphson
1. Implement the Newton-Raphson algorithm
2. Test with $g(x)$.

   $$
   g(x) = 2 \cos x - x
   $$
   
3. Bonus: test performance of `newton_raphson()` against `bisection()`.


```python
def g(x):
    return 2*np.cos(x) - x
```


```python
xvals = np.linspace(0, 7, 30)
plt.plot(xvals, np.zeros_like(xvals), 'k--')
plt.plot(xvals, g(xvals))
```


```python
def newton_raphson(f, x, h=3e-1, Nmax=100, eps=1e-14):
    for iteration in range(Nmax):
        fx = f(x)
        if np.abs(fx) < eps:
            break
        df = (f(x + h/2) - f(x - h/2))/h
        Delta_x = -fx/df
        x += Delta_x
    else:
        print("Newton-Raphson: no root found after {0} iterations (eps={1}); "
              "best guess is {2} with error {3}".format(Nmax, eps, x, fx))
        x = None
    return x
```


```python
newton_raphson(g, 0)
```




    1.0298665293222589




```python

```

ca7df84862923f8c3c48343b9435bc090404f99c
ipynb
Jupyter Notebook
13_root_finding/.ipynb_checkpoints/13-Root-finding-checkpoint.ipynb
nachrisman/PHY494
bac0dd5a7fe6f59f9e2ccaee56ebafcb7d97e2e7
CC-BY-4.0
13_root_finding/.ipynb_checkpoints/13-Root-finding-checkpoint.ipynb
nachrisman/PHY494
bac0dd5a7fe6f59f9e2ccaee56ebafcb7d97e2e7
CC-BY-4.0
13_root_finding/.ipynb_checkpoints/13-Root-finding-checkpoint.ipynb
nachrisman/PHY494
bac0dd5a7fe6f59f9e2ccaee56ebafcb7d97e2e7
CC-BY-4.0
Qwen/Qwen-72B
1. YES
2. YES