<a href="https://colab.research.google.com/github/stephenbeckr/numerical-analysis-class/blob/master/Demos/Ch2_Intersection_GraphingCalculator.ipynb" target="_parent"></a>

# Finding the intersection of two functions

This is the **master** branch version. We suggest you first try the [**student**](https://github.com/stephenbeckr/numerical-analysis-class/blob/student/Demos/Ch2_Intersection_GraphingCalculator.ipynb) branch version, which has open-ended questions.

By the way, if the "Open in Colab" link doesn't work, try right-clicking the link and choose "open in new tab".


### Finding intersections of two functions
This is something you might have done in high school on a graphing calculator (the particular functions we're using are inspired by this [post about using the TI-84+](https://www.dummies.com/education/graphing-calculators/how-to-find-points-of-intersection-on-the-ti-84-plus/)).



Let's find the intersection of
$$ f(x) = .2(x-2)x(x+4)  \quad\text{and}\quad g(x) = .5 x $$

First, let's plot the functions:
(for plotting, here's a list of [style sheets](https://matplotlib.org/3.3.1/gallery/style_sheets/style_sheets_reference.html) to use with `plt.style.use`, and [general matplotlib cheatsheets](https://github.com/matplotlib/cheatsheets); if you want to make the axes through the origin, like in the TI-84 picture, this [stackoverflow post](https://stackoverflow.com/a/25689340) gives several methods)


```
import numpy as np


# basic plots
%matplotlib inline  

# allows interative plots, so you can zoom in!
# However, can be a bit buggy
#%matplotlib notebook

#%matplotlib ipympl # see https://medium.com/@1522933668924/using-matplotlib-in-jupyter-notebooks-comparing-methods-and-some-tips-python-c38e85b40ba1

import matplotlib.pyplot as plt


plt.style.use('seaborn-notebook') # just a personal preference for the style

# == First, define the function f
f = lambda x : .2*(x-2)*x*(x+4)


# == Now, define g
g = lambda x : .5*x


# == Now, plot

x   = np.linspace(-10,10,100)
```


```
fig, ax = plt.subplots()
ax.plot( x, f(x) )
ax.plot( x, g(x) )
ax.set_ylim(-10,10)
plt.xticks(range(-10,11,5))
plt.yticks(range(-10,11,5))

# Adjust it to look a bit more like the TI-84 image
# Alternatively, see https://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.axvline
ax.spines['left'].set_position('zero')
ax.spines['right'].set_color('none')
ax.spines['bottom'].set_position('zero')
ax.spines['top'].set_color('none')
ax.minorticks_on()
fig.show() # Optional in jupyter
```

## Now, let's find all the roots (it seems we have three of them, looking at the plots).

Use the builtin root finder from `scipy` (see [`scipy.optimize` documentation](https://docs.scipy.org/doc/scipy/reference/optimize.html#root-finding)) called `root_scalar`

You can either specify the method you want, or it will pick a default method for you, based on what you supply, e.g.,
- if you supply an **interval** containing the root, using `bracket = [a,b]`, then it will use one of Brent's methods; you can also request using the bisection method
- if you supply an **initial guess** using `x0 = 1.2` (for example), then it will try something like Newton's method. However, you need to give it a bit more information:
  - to actually use **Newton's method** you have to specify a function that computes the derivative, using `fprime = [insert name-of-derivative-function]`
  - you can use the **secant method** which doesn't need the derivative, but then you have to specify a second initial point via `x1 = 1.3` (for example)


```
from scipy.optimize import root_scalar

output = root_scalar( lambda x : f(x) - g(x), x0 = -4, x1 = -4.1)
r1     = output.root
print('Found first root at x={} and y={}\n'.format(r1,f(r1)) ) # see https://pyformat.info/ for print statement formatting
output
```

    Found first root at x=-4.3911649915626345 and y=-2.19558249578132
    





          converged: True
               flag: 'converged'
     function_calls: 7
         iterations: 6
               root: -4.3911649915626345




```
output = root_scalar( lambda x : f(x) - g(x), x0 = 0, x1 = .1)
r2     = output.root
print('Found second root at x={} and y={}\n'.format(r2,f(r2)) ) # see https://pyformat.info/ for print statement formatting
output
```

    Found second root at x=0.0 and y=-0.0
    





          converged: True
               flag: 'converged'
     function_calls: 3
         iterations: 2
               root: 0.0




```
output = root_scalar( lambda x : f(x) - g(x), x0 = 2.1, x1 = 2.2)
r3     = output.root
print('Found third root at x={} and y={}\n'.format(r3,f(r3)) ) # see https://pyformat.info/ for print statement formatting
output
```

    Found third root at x=2.3911649915626345 and y=1.1955824957813186
    





          converged: True
               flag: 'converged'
     function_calls: 7
         iterations: 6
               root: 2.3911649915626345




```
ax.plot( r1, f(r1), 'o',MarkerSize=12 )
ax.plot( r2, f(r2), 'o',MarkerSize=12 )
ax.plot( r3, f(r3), 'o',MarkerSize=12 )
fig
```

### Compare accuracy with `roots`

Be careful: despite a similar sounding name, [`numpy.roots`](https://numpy.org/doc/stable/reference/generated/numpy.roots.html) is **only for polynomials**.

With a bit of pen-and-paper manipulation, you can cast the above intersection problem as one of finding the roots of a polynomial.  Do this, and then find the roots via `numpy.roots`.

- How can you evaluate the accuracy of your answer?
- Which method (`scipy.optimize.root_scalar` vs `numpy.roots`) is more accurate?



```
# Make the polynomial we want to find roots of (i.e., find the coefficients)

coeff = np.poly( (2,0,-4) )  # coefficients of (x-2)*x*(x+4)
# and don't forget to multiply by 0.2
coeff = 0.2*coeff
# and subtract the g(x) function, so -.5*x.  This is the 2nd-to-last coefficient
coeff[-2] -=.5
coeff
```




    array([ 0.2,  0.4, -2.1,  0. ])




```
r_numpy = np.roots( coeff )
r_numpy
```




    array([-4.39116499,  2.39116499,  0.        ])




```
r_scipy = np.array( (r1,r3,r2 ))  # changed ordering to match np.roots output

residual = lambda x : abs( f(x) - g(x) )
for r_n, r_s in zip(r_numpy,r_scipy):
    print("At root {:5.2f}, the residual using np.roots is {:.2e}, and using root_scalar is {:.2e}".format( r_n, residual(r_n), residual(r_s)) )
```

    At root -4.39, the residual using np.roots is 2.66e-15, and using root_scalar is 2.66e-15
    At root  2.39, the residual using np.roots is 2.22e-16, and using root_scalar is 1.33e-15
    At root  0.00, the residual using np.roots is 0.00e+00, and using root_scalar is 0.00e+00


Conclusion: in this example, both root finding libraries work just fine

## Let's repeat, for a different function

Again, let's pick a polynomial since we can find the true roots by hand

Let's do $$ f(x) = (x-1) \cdot (x-(1+10^{-8}))^2. $$

Do we expect any conditioning issues?

Yes!  This is almost a double-root at $x=1$, so very shallow slope.  If we perturb $f$ upward a bit, it won't have any roots at all!

### Try to find the roots using the polynomial root-finding routines (`np.roots`)


```

trueRoots = (1,1+1e-8)
coeff = np.poly( trueRoots  )
f     = lambda x : np.polyval( coeff, x )

r_numpy = np.roots( coeff )
r_numpy
```




    array([1.+1.49011612e-08j, 1.-1.49011612e-08j])



How did it go?

We got imaginary numbers! This means there were no real roots. But there are! It's just a hard problem, and very ill-conditioned.

The plot of the function will help us


```
fig, ax = plt.subplots(1,2);
x = np.linspace(0,2,100)
ax[0].plot( x, f(x) );

# Zoom in a bit
x = np.linspace(1-1e-8,1+5e-8,100);
ax[1].plot( x, f(x) );
ax[1].spines['bottom'].set_position('zero')
ax[1].spines['top'].set_color('none')
ax[1].spines['right'].set_color('none')
f( x[9] )

# from sympy import init_printing
# init_printing()
# list(zip( range(20), f( x[:20] ).tolist() ))

for xx, fx in zip( x[:20], f( x[:20] ) ):
  print("At x = 1 + {:5.2e}, f(x) is {:5.2e}".format(xx-1,fx) )
```

The table of values shows us the problem!  Within our limited precision, **many** values of x give us a root, since the nearest floating point number is $0$

## Let's try finding the roots via a generic rootfinding method (i.e., ignoring the fact that we have a polynomial)

Call [`root_scalar`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.root_scalar.html#scipy.optimize.root_scalar). 


```
coeff_deriv = np.polyder( coeff )
deriv   = lambda x : np.polyval( coeff_deriv, x )
root_scalar( f, x0 = 1.1, fprime = deriv )
```




          converged: True
               flag: 'converged'
     function_calls: 46
         iterations: 23
               root: 1.0000000162008909



Let's repeat, but with a different evaluation of $f$ (taking advantage that we already know it in factored form)


```
f       = lambda x : (x-trueRoots[0])*(x-trueRoots[1]) # different
coeff_deriv = np.polyder( coeff )
deriv   = lambda x : np.polyval( coeff_deriv, x )
out = root_scalar( f, x0 = 1.1, fprime = deriv )
out
```




          converged: True
               flag: 'converged'
     function_calls: 46
         iterations: 23
               root: 1.0000000176119155




```
# What is the relative error?
abs( out.root - trueRoots[1] )/abs(trueRoots[1])
```


```
# We have about 8 digits of accuracy, roughly
```

fe009c6e28f9575bda59dc17631f97c447a30601
ipynb
Jupyter Notebook
Demos/Ch2_Intersection_GraphingCalculator.ipynb
skhadem/numerical-analysis-class
a022fc2e73254800a1e193f94280446223ef01ae
BSD-3-Clause
2020-08-25T19:11:46.000Z
2021-11-12T21:49:44.000Z
Demos/Ch2_Intersection_GraphingCalculator.ipynb
sabrinazhengliu/numerical-analysis-class
1585b68d3f3c375259507bf18f03ba0332851edb
BSD-3-Clause
2020-09-01T21:44:12.000Z
2020-09-01T21:44:12.000Z
Demos/Ch2_Intersection_GraphingCalculator.ipynb
sabrinazhengliu/numerical-analysis-class
1585b68d3f3c375259507bf18f03ba0332851edb
BSD-3-Clause
2020-08-25T21:25:17.000Z
2021-03-03T04:14:26.000Z
Qwen/Qwen-72B
1. YES
2. YES