** Counting roots
:PROPERTIES:
:categories: Nonlinear algebra
:date: 2013/02/27 10:13:59
:updated: 2013/02/27 14:27:48
:END:
[[http://matlab.cheme.cmu.edu/2011/09/10/counting-roots/][Matlab post]]
The goal here is to determine how many roots there are in a nonlinear function we are interested in solving. For this example, we use a cubic polynomial because we know there are three roots.
$$f(x) = x^3 + 6x^2 - 4x -24$$
*** Use roots for this polynomial
This ony works for a polynomial, it does not work for any other nonlinear function.
#+BEGIN_SRC python
import numpy as np
print np.roots([1, 6, -4, -24])
#+END_SRC
#+RESULTS:
: [-6. 2. -2.]
Let us plot the function to see where the roots are.
#+BEGIN_SRC python
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(-8, 4)
y = x**3 + 6 * x**2 - 4*x - 24
plt.plot(x, y)
plt.savefig('images/count-roots-1.png')
#+END_SRC
#+RESULTS:
[[./images/count-roots-1.png]]
Now we consider several approaches to counting the number of roots in this interval. Visually it is pretty easy, you just look for where the function crosses zero. Computationally, it is tricker.
*** method 1
Count the number of times the sign changes in the interval. What we have to do is multiply neighboring elements together, and look for negative values. That indicates a sign change. For example the product of two positive or negative numbers is a positive number. You only get a negative number from the product of a positive and negative number, which means the sign changed.
#+BEGIN_SRC python
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(-8, 4)
y = x**3 + 6 * x**2 - 4*x - 24
print np.sum(y[0:-2] * y[1:-1] < 0)
#+END_SRC
#+RESULTS:
: 3
This method gives us the number of roots, but not where the roots are.
*** Method 2
Using events in an ODE solver python can identify events in the solution to an ODE, for example, when a function has a certain value, e.g. f(x) = 0. We can take advantage of this to find the roots and number of roots in this case. We take the derivative of our function, and integrate it from an initial starting point, and define an event function that counts zeros.
$$f'(x) = 3x^2 + 12x - 4$$
with f(-8) = -120
#+BEGIN_SRC python
import numpy as np
from pycse import odelay
def fprime(f, x):
return 3.0 * x**2 + 12.0*x - 4.0
def event(f, x):
value = f # we want f = 0
isterminal = False
direction = 0
return value, isterminal, direction
xspan = np.linspace(-8, 4)
f0 = -120
X, F, TE, YE, IE = odelay(fprime, f0, xspan, events=[event])
for te, ye in zip(TE, YE):
print 'root found at x = {0: 1.3f}, f={1: 1.3f}'.format(te, ye)
#+END_SRC
#+RESULTS:
: root found at x = -6.000, f=-0.000
: root found at x = -2.000, f=-0.000
: root found at x = 2.000, f= 0.000