Course Webpage

Introduction to python for numerical computing II

Horner algorithm

We import matplotlib.pyplot and numpy modules and the functions of numpy for polynomial.

In [3]:
import numpy as np
import matplotlib.pyplot as plt
import numpy.polynomial.polynomial as pol

Given a polynomial

$$ P(x) = p_0+p_1\,x+p_2\,x^2+p_3\,x^3+p_4\,x^4+p_5\,x^5 $$

For the next example, we will use as as coefficients of this polynomial

In [4]:
p = np.array([1., -1, 2, -3,  5, -2])
p0, p1, p2, p3, p4, p5 = p

So the polynomial is

$$ P(x) = 1-x+2\,x^2-3\,x^3+5\,x^4-2\,x^5 $$

We can define this polynomial and its derivatives as lambda functions

In [6]:
g = lambda x: p0 + p1*x + p2*x**2 + p3*x**3 + p4*x**4 + p5*x**5
d1g = lambda x: p1 + 2*p2*x + 3*p3*x**2 + 4*p4*x**3 + 5*p5*x**4
d2g = lambda x: 2*p2 + 6*p3*x + 12*p4*x**2 + 20*p5*x**3

We can plot them

In [7]:
a = 0.; b = 1.
x = np.linspace(a,b)
plt.plot(x,0*x,'k')
plt.plot(x,g(x), label = 'g')
plt.plot(x,d1g(x), label = 'd1g')
plt.plot(x,d2g(x), label = 'd2g')
plt.legend()
plt.show()

There is a better way of working with polynomials: once we have stored the coefficients in a numpy array we do not need to use lambda functions.

We can compute the polynomial value at any point (or many points simultaneously) using the numpy command polyval.

In [9]:
x0 = -0.5
print('P value at point       ', x0)
print('With polyval:          ', pol.polyval(x0, p))
print('With lambda function g:', g(x0))

P value at point        -0.5
With polyval:           2.75
With lambda function g: 2.75

Let us implement the algorithm used by polyval, the Horner algorithm, for unidimensional arrays. First, let us calculate how many additions and multiplications we would need to compute the polynomial at a point $x_0$.

$$ P(x_0) = p_0+p_1\,x_0+p_2\,x_0^2+p_3\,x_0^3+p_4\,x_0^4+p_5\,x_0^5 $$

We need 1 multiplication for the second term, 2 for the third one, 3 for the fourth one, 4 for the fith one and 5 for the last one. Also, we must perform 5 additions. That is,

  • $1+2+3+4+5=15$ multiplications.
  • $5$ additions.

Let us see how we compute the polynomial with Horner algorithm.

If we apply Ruffini to the polynomial with $x_0 = 1$

$$ \begin{array}{c|rrrrrrl} & -2 & 5 & -3 & 2 & -1 & 1 & \\ 1 & & -2 & 3 & 0 & 2 & 1 & \\ \hline & -2 & 3 & 0 & 2 & 1 & \fbox{2}&=P(1) \end{array} $$

If $$P(x)= 1-x+2x^2-3x^3+5x^4-2x^5$$ and $$Q(x)=1+2x+3x^3-2x^4,$$ then $$P(x) = Q(x)(x-1)+2$$ and it follows $$P(1)= Q(1)(1-1)+2,$$ that is, $$P(1)=2.$$

And the operations performed are

$$ \begin{array}{c|rrrrrrl} & -2 & 5 & -3 & 2 & -1 & 1 & \\ 1 & & -2 & 3 & 0 & 2 & 1 & \\ \hline & -2 & 3 & 0 & 2 & 1 & \fbox{2}&=P(1) \end{array} $$

That is

$$ \begin{array}{c|ccccccl} & p_5 & p_4 & p_3 & p_2 & p_1 & p_0 & \\ x_0 & & q_5\,x_0 & q_4\,x_0 & q_3\,x_0 & q_2\,x_0 & q_1\,x_0 & \\ \hline & q_5 & q_4 & q_3 & q_2 & q_1 & \fbox{$q_0$}&=P(x_0) \end{array} $$$$ \begin{array}{l} q_5 = p_5\\ q_4 = p_4 + q_5\,x_0\\ q_3 = p_3 + q_4\,x_0\\ q_2 = p_2 + q_3\,x_0\\ q_1 = p_1 + q_2\,x_0\\ q_0 = p_0 + q_1\,x_0\\ \end{array} $$

That is, $5$ multiplications and $5$ additions. In general, for an $n$ degree polynomial, $n$ multiplications and $n$ additions will be performed.

Exercises

Exercise 1

Write a function that implements the Horner algorithm. The input arguments will be an unidimensional numpy array (vector) with the $P$ polynomial coefficients and the point $x_0$. That is, horner(x0, p). The output arguments will be a unidimensional numpy array that contains the coefficients of $Q$ and we will call quotient, and a float variable that will be $P\,(x_0)$ and we will call it remainder. Print the coefficients of $Q$ and $P\,(x_0)$. Check that the value of $P\,(x_0)$ is correct using pol.polyval(x0,p)

Check the results for the polynomial $P_0$ and the point $x_0$

p0 = np.array([1.,2,1])
x0 = 1.

Also, check the results for the polynomial $P_1$ and the point $x_1$

p1 = np.array([1., -1, 2, -3,  5, -2])
x1 = 1.

Also, check the results for the polynomial $P_2$ and the point $x_2$

p2 = np.array([1., -1, -1, 1, -1, 0, -1, 1])
x2 = -1.

Hints:

  • The polynomial $P(x) = p_0 + p_1 x + p_2 x^2 + p_3 x^3$ is stored p = np.array([p0,p1,p2,p3])
  • The number of elements of a vector p can be obtained with n = len(p).
  • We can initialize a vector q of the same size than p with q = np.zeros_like(p) or q = np.zeros(n).
  • range(n1,n2,-1) creates a consecutive integer sequence from n1 to n2 (without reaching n2) with step -1.
  • The last element of p is es p[n-1] or p[-1].
  • If we want to extract all the elements of a p but the first we can write v[1:].
  • The last line of the function could be
    return quotient, remainder
  • And we can call the function with
    q, r = horner(x0, p)

To make easier the importation of the horner function in the exercise 3, the structure of the code of this exercise will be

import numpy as np
import numpy.polynomial.polynomial as pol 

def horner(x0,p):
    ...  # <--- write here your code
    return quotient, remainder

def main():
    p0 = np.array([1.,2,1])
    x0 = 1.
    q, r = horner(x0,p0)
    rp   = pol.polyval(x0,p0) 

    print('\nQ coefficients = ', q)
    print('P0(1)        = ', r)
    print('With polyval = ', rp)

    ... # <- examples with polynomials p1 and p2

if __name__ == "__main__":
    main()
In [11]:
%run Exercise1.py

Q coefficients =  [3. 1.]
P0(1)        =  4.0
With polyval =  4.0


Q coefficients =  [ 1.  2.  0.  3. -2.]
P1(1)        =  2.0
With polyval =  2.0


Q coefficients =  [ 4. -5.  4. -3.  2. -2.  1.]
P2(-1)       =  -3.0
With polyval =  -3.0

Exercise 2

Modify the previous function in such a way that now, the input argument is a vector x = np.linspace(-1,1), instead of a point x0.

The function will be HornerV(x,p). Now, the output argument will be a vector y that contains the values of polynomial $P$ for the points in x . Use an external loop that repeats the process for each point in x and that stores every result in y. Check that it works plotting the polynomial in the interval $[-1,1].$

Plot the polynomial $P$ of coefficients

p = np.array([1., -1, 2, -3, 5, -2])

Plot also the polynomial $R$ of coefficients

r = np.array([1., -1, -1, 1, -1, 0, -1, 1])

Hint:

  • We can initialize a vector of zeros y with the same size as x with y = np.zeros_like(x).
In [13]:
%run Exercise2.py

We can also do it with polyval command.

In [17]:
x = np.linspace(-1,1)

plt.figure()
plt.plot(x,pol.polyval(x,r))
plt.plot(x,0*x,'k')
plt.title('R')
plt.show()

Taylor series for a polynomial of degree 5

\begin{gather*} P(x) = p_0+p_1x+p_2x^2+p_3x^3+p_4x^4+p_5x^5 \\[20pt] P(x) = P(x_0) + \frac{P'(x_0)}{1!}(x-x_0) + \frac{P''(x_0)}{2!}(x-x_0)^{2} + \frac{P'''(x_0)}{3!}(x-x_0)^{3} + \frac{P^{(4)}(x_0)}{4!}(x-x_0)^{4} + \frac{P^{(5)}(x_0)}{5!}(x-x_0)^{5} \\[20pt] P(x) = P(x_0) + (x-x_0) \underset{Q_1(x)}{\left[ \underbrace{\frac{P'(x_0)}{1!}+(x-x_0) \underset{Q_2(x)}{\left[ \underbrace{\frac{P''(x_0)}{2!}+(x-x_0) \underset{Q_3(x)}{\left[ \underbrace{\frac{P'''(x_0)}{3!}+(x-x_0) \underset{Q_4(x)}{\left[ \underbrace{\frac{P^{(4)}(x_0)}{4!}+(x-x_0) \underset{Q_5(x)}{\left[ \underbrace{\frac{P^{(5)}(x_0)}{4!} }\right]} }\right]} }\right]} }\right]} }\right]} \end{gather*}

thus, if we apply the Horner algorithm to the resultant polynomials, we can obtain the higher order derivative values at the point. For instance, for $x_0=1$ and $P(x) = 1-x+2\,x^2-3\,x^3+5\,x^4-2\,x^5$

$$ \begin{array}{c|rrrrrrl} & -2 & 5 & -3 & 2 & -1 & 1 &\\ 1 & & -2 & 3 & 0 & 2 & 1 & \\ \hline & -2 & 3 & 0 & 2 & 1 & \fbox{2}&=P(1)\,/\,0!\\ 1 & & -2 & 1 & 1 & 3 & & \\ \hline & -2 & 1 & 1 & 3 & \fbox{4} & &=P\,'(1)\,/\,1! \\ 1 & & -2 & -1 & 0 & & & \\ \hline & -2 & -1 & 0 & \fbox{3} & & &=P\,''(1)\,/\,2! \\ 1 & & -2 & -3 & & & & \\ \hline & -2 & -3 & \fbox{$-3$} & & & &=P\,'''(1)\,/\,3! \\ 1 & & -2 & & & & & \\ \hline & -2 & \fbox{$-5$} & & & & &=P\,^{(4)}(1)\,/\,4! \\ 1 & & & & & & & \\ \hline & \fbox{$-2$} & & & & & &=P\,^{(5)}(1)\,/\,5! \\ \end{array} $$

Exercise 3

Using the function from exercise 1, compute and print the higher order derivatives of $P(x)$ at $x_0=1$ using a function polDer(x0,p). Do not use a factorial function: generate the factorial value similarly as in the previous lab for McLaurin series.

Check the results for the polynomial $P$ and the point $x_0$

(a) Write the remainders of dividing $P$ again and again by $(x-x_0).$ Same for $R$ and $x_1$.

p = np.array([1., -1, 2, -3,  5, -2])
x0 = 1.

r = np.array([1., -1, -1, 1, -1, 0, -1, 1])
x1 = -1.

(b)

Write the successive derivatives of $P$ at the point $x_0.$ Same for $R$ and $x_1$. To obtain them we will multiply each remainder byt the corresponding factorial because

$$ r_i=\frac{P\,^{(i)}(x_0)}{i!} \quad \Longrightarrow \quad P\,^{(i)}(x_0) = r_i \times i!$$

Hints:

  • If we want to duplicate a numpy array, write v2 = np.copy(v1) because if we code v2 = v1 the two arrays are stored in the same location. Then, if we modify one element in one vector, this element is also changed in the other one.
  • You can use np.set_printoptions(suppress = True) to avoid exponential notation. Write it after importing numpy.
  • If we want to import the horner function from exercise 1, we can write at the top of the file
    from Exercise1 import horner
In [20]:
%run Exercise3a.py

Remainders of dividing P again and again by (x-x0)
[ 2.  4.  3. -3. -5. -2.]


Remainders of dividing R again and again by (x-x1)
[ -3.  21. -46.  60. -51.  27.  -8.   1.]
In [22]:
%run Exercise3b.py

Derivatives of P in x0 = 1
[   2.    4.    6.  -18. -120. -240.]


Derivatives of R in x1 = -1
[   -3.    21.   -92.   360. -1224.  3240. -5760.  5040.]

Proposed exercises

Let be a polynomial with $a_n=1$ and the independent term nonzero and where all its roots are integers and single.

If we apply Horner's algorithm once and again and we try different divisors of the independent term, we can find all the roots of the polynomial.

For example, the independent term of the polynomial

$$ P(x) = 18+9x-20x^2-10x^3+2x^4+x^5 $$

is $18$. Its divisors are $\{1,2,3,6,9,18\}$ and possible integer roots are $\{1,-1,2,-2,3,-3,6,-6,9,-9,18,-18\}.$

If we apply horner

$$ \begin{array}{c|rrrrrrl} & 1 & 2 & -10 & -20 & 9 & 18 &\\ 1 & & 1 & 3 & -7 & -27 & -18 & \\ \hline & 1 & 3 & -7 & -27 & -18 & \fbox{0}\\ -1 & & -1 & -2 & 9 & 18 & & \\ \hline & 1 & 2 & -9 & -18 & \fbox{0} & & \\ -2 & & -2 & 0 & 18 & & & \\ \hline & 1 & 0 & -9 & \fbox{0} & & & \\ 3 & & 3 & 9 & & & & \\ \hline & 1 & 3 & \fbox{0} & & & & \\ -3 & & -3 & & & & & \\ \hline & 1 & \fbox{0} & & & & &\\ \end{array} $$

we can conclude that

$$ P(x) = (x-1)(x+1)(x+2)(x-3)(x+3) $$

Exercise 4

(a)

Write a function divisors(m) that computes the integer divisors of a given integer number m. The output argument will be a numpy array that contains its divisors with positive and negative sign. For example, for m = 18 the output will be the numpy array div = [1,-1,2,-2,3,-3,6,-6,9,-9,18,-18].

Hints:

  • The indices must be integer. Also, we want a positive m . So we initialize the vector that will contain the divisors
    m = abs(int(m))
div = np.zeros(2*m)
    and after storing the divisors, we chop to the appropiate size n, for example, with div = div[:n].
  • To compute the remainder of the integer division, np.remainder(x1,x2) can be used. It is equivalent to the Python modulus operator x1 % x2.

(b)

Write a function roots(p) which computes the single integer roots of a polynomial and returns them in a numpy array. Use the functions divisors and horner.

Taking into account that all the roots of the following polynomials are integers and simple, compute the roots of the polynomials

p0 = np.array([-1., 0, 1])
p1 = np.array([8., -6, -3, 1])
p2 = np.array([15., -2, -16, 2, 1])
p3 = np.array([60., 53, -13, -5, 1])  
p4 = np.array([490., 343, -206, -56, 4, 1])

Hints:

  • Take into account that the number of roots of a polynomial is the degree of the polynomial.
  • We can proceed the following way: once we have found a root (the remainder must be zero)
    1. We store the x0 (a divisor of the constant term of the polynomial) as a root.
    2. For the next step, we use the last quotient as the next polynomial to divide by $(x-x_0)$.
    3. We try with the next divisor of the constant term of the polynomial.
In [1]:
%run Exercise4.py

(a)


Divisors of 6
[ 1. -1.  2. -2.  3. -3.  6. -6.]

Divisors of 18
[  1.  -1.   2.  -2.   3.  -3.   6.  -6.   9.  -9.  18. -18.]

Divisors of 20
[  1.  -1.   2.  -2.   4.  -4.   5.  -5.  10. -10.  20. -20.]


(b)


Roots of p0

[ 1. -1.]

Roots of  p1

[ 1. -2.  4.]

Roots of p2

[ 1. -1.  3. -5.]

Roots of p3

[-1. -3.  4.  5.]

Roots of p4

[-1.  2. -5.  7. -7.]

Given the polynomial with $a_n=1$

$$ P(x) = 8-22x+17x^2+x^3-5x^4+x^5 $$

$a_0=8$ and its divisors are $\{1,2,4,8\}$ and possible integer roots are $\{1,-1,2,-2,4,-4,8,-8\}$

If we apply horner

$$ \begin{array}{c|rrrrrrl} & 1 & -5 & 1 & 17 & -22 & 8 &\\ 1 & & 1 &-4 & -3 & 14 & -8 & \\ \hline & 1 & -4 & -3 & 14 & -8 & \fbox{0}\\ 1 & & 1 & -3 & -6 & 8 & & \\ \hline & 1 & -3 & -6 & 8 & \fbox{0} & & \\ 1 & & 1 & -2 & -8 & & & \\ \hline & 1 & -2 & -8 & \fbox{0} & & & \\ -2 & & -2 & 8 & & & & \\ \hline & 1 & -4 & \fbox{0} & & & & \\ 4 & & 4 & & & & & \\ \hline & 1 & \fbox{0} & & & & &\\ \end{array} $$

we can conclude that

$$ P(x) = (x-1)(x-1)(x-1)(x+2)(x-4) $$

that is

$$ P(x) = (x-1)^3(x+2)(x-4) $$

Exercise 5

Modify the Exercise 4 in such a way that we can compute integer multiple roots of a polynomial. If a root is multiple, it will appear in the array as many times as its multiplicity.

Compute the integer roots of

p1 = np.array([8., -22, 17, 1, -5, 1])
p2 = np.array([-135., 378, -369, 140, -9, -6, 1])
p3 = np.array([96., 320, 366, 135, -30, -24, 0, 1])
p4 = np.array([280., 156, -350, -59, 148, -26, -6, 1])
In [2]:
%run Exercise5.py

Roots of p1

[ 1.  1.  1. -2.  4.]

Roots of p2

[ 1.  1.  3.  3.  3. -5.]

Roots of p3

[-1. -1. -1. -2. -3.  4.  4.]

Roots of p4

[-1. -1.  2.  2.  2. -5.  7.]

If we apply Horner

$$ \begin{array}{c|ccccccl} & p_5 & p_4 & p_3 & p_2 & p_1 & p_0 & \\ x_0 & & q_5\,x_0 & q_4\,x_0 & q_3\,x_0 & q_2\,x_0 & q_1\,x_0 & \\ \hline & q_5 & q_4 & q_3 & q_2 & q_1 & \fbox{$q_0$}&=P(x_0) \end{array} $$$$ \begin{array}{l} q_5 = p_5\\ q_4 = p_4 + q_5\,x_0\\ q_3 = p_3 + q_4\,x_0\\ q_2 = p_2 + q_3\,x_0\\ q_1 = p_1 + q_2\,x_0\\ q_0 = p_0 + q_1\,x_0\\ \end{array} $$

If we do not need to keep the intermediate values, only the last one, we can use a float variable, instead of a numpy array

$$ \begin{array}{c|ccccccl} & p_5 & p_4 & p_3 & p_2 & p_1 & p_0 & \\ x_0 & & y_0\,x_0 & y_0\,x_0 & y_0\,x_0 & y_0\,x_0 & y_0\,x_0 & \\ \hline & y_0 & y_0 & y_0 & y_0 & y_0 & \fbox{$y_0$}&=P(x_0) \end{array} $$$$ \begin{array}{l} y_0 = p_5\\ y_0 = p_4 + y_0\,x_0\\ y_0 = p_3 + y_0\,x_0\\ y_0 = p_2 + y_0\,x_0\\ y_0 = p_1 + y_0\,x_0\\ y_0 = p_0 + y_0\,x_0\\ \end{array} $$

And we have the advantage that it works both for float values as $x_0,$ $y_0$ and numpy arrays $x,$ $y.$ That is, we can use

$$ \begin{array}{l} y = p_5\\ y = p_4 + y\,x\\ y = p_3 + y\,x\\ y = p_2 + y\,x\\ y = p_1 + y\,x\\ y = p_0 + y\,x\\ \end{array} $$

Exercise 6

Repeat Exercise 2 vectorizing the process as follows:

  • We do not add any loop to the function horner(x0,p). The function will be now hornerVect(x,p).
  • We will use a vector y that, in the last step, will contain all the values of the polynomial for the points in the vector x.
  • The value of y in the first step is the same for all the x points, p[n-1] with n = len(p), so we can use y = p[n-1] or y = p[-1].

Check, without plotting, with the polynomial pO and the points x0 and print y

p0 = np.array([1.,2,1])
x0 = np.array([1.,-1])

Plot in $[-1,1]$ the polynomials

p = np.array([1., -1, 2, -3, 5, -2])
r = np.array([1., -1, -1, 1, -1, 0, -1, 1])

Using the vector x = np.linspace(-1,1,1000000) compare the execution time with this code and that from Exercise 2.

Check the results for the polynomial $P$

p = np.array([1., -1., 2., -3., 5., -2.])
In [7]:
%run Exercise6.py

y =  [4. 0.]



Time without vectorization =  1.1784281730651855
Time with vectorization    =  0.009869098663330078

Exercise 7

Modify the script in Exercise 3:

  • Now, the input argument must be a vector x = np.linspace(0,1) instead of a value x0. The function will be derivativesH(x,p).
  • The output argument will be an Y matrix that will contain, in the first column, the values of $P$ at the points in x, in the second column, the values of $P'$ at the same points, and so on until $n$-th derivative, being $n$ the $P$ degree.
  • Using this Y matrix, plot the polynomial and its first and second derivatives in $[-1,1]$.

Hints:

  • Matrix Y can be initialized Y = np.zeros((m,n)) where m = len(x) and n = len(p).
  • If we want to extract a colum k from a matrix Y we can write Y[:,k].

Plot $P$, $P'$ and $P''$

p = np.array([1., -1., 2., -3.,  5., -2.])
In [8]:
%run Exercise7.py