5

it is commonly an easy task to build an n-th order polynomial and find the roots with numpy:

import numpy
f = numpy.poly1d([1,2,3])
print numpy.roots(f)
array([-1.+1.41421356j, -1.-1.41421356j])

However, suppose you want a polynomial of type:

f(x) = a*(x-x0)**0 + b(x-x0)**1 + ... + n(x-x0)**n

Is there a simple way to construct a numpy.poly1d type function and find the roots ? I've tried scipy.fsolve but it is very unstable as it depends highly on the choice of the starting values in my particular case.

Thanks in advance Best Regards rrrak

EDIT: Changed "polygon"(wrong) to "polynomial"(correct)

rrrak
  • 3
  • 1
rrrak
  • 51
  • 1
  • 2

2 Answers2

4

First of all, surely you mean polynomial, not polygon?

In terms of providing an answer, are you using the same value of "x0" in all the terms? If so, let y = x - x0, solve for y and get x using x = y + x0.

You could even wrap it in a lambda function if you want. Say, you want to represent

f(x) = 1 + 3(x-1) + (x-1)**2

Then, 

>>> g = numpy.poly1d([1,3,1])
>>> f = lambda x:g(x-1)
>>> f(0.0)
-1.0

The roots of f are given by:

f.roots = numpy.roots(g) + 1
Pascal Bugnion
  • 4,878
  • 1
  • 24
  • 29
3

In case x0 are different by power, such as:

f(x) = 3*(x-0)**0 + 2*(x-2)**1 + 3*(x-1)**2 + 2*(x-2)**3

You can use polynomial operation to calculate the finally expanded polynomial:

import numpy as np
import operator

ks = [3,2,3,2]
offsets = [0,2,1,2]

p = reduce(operator.add, [np.poly1d([1, -x0])**i * c for i, (c, x0) in enumerate(zip(ks, offsets))])

print p

The result is:

   3     2
2 x - 9 x + 20 x - 14
HYRY
  • 94,853
  • 25
  • 187
  • 187
  • this sounds a good solution to reduce polynomials but why not simple: `sum([np.poly1d([1, -x0])**i * c for i, (c, x0) in enumerate(zip(ks, offsets))])` – dashesy Dec 26 '15 at 02:28