Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

numpy.poly1d , root-finding optimization, shifting polynom on x-axis

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)

like image 1000
rrrak Avatar asked Jul 22 '11 14:07

rrrak


2 Answers

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
like image 110
Pascal Bugnion Avatar answered Sep 17 '22 19:09

Pascal Bugnion


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
like image 40
HYRY Avatar answered Sep 17 '22 19:09

HYRY