** Polynomials in python
:PROPERTIES:
:categories: math, polynomials
:date: 2013/01/22 09:00:00
:updated: 2013/02/27 14:53:59
:END:
[[http://matlab.cheme.cmu.edu/2011/08/01/polynomials-in-matlab/][Matlab post]]
Polynomials can be represented as a list of coefficients. For example, the polynomial $4*x^3 + 3*x^2 -2*x + 10 = 0$ can be represented as [4, 3, -2, 10]. Here are some ways to create a polynomial object, and evaluate it.
#+BEGIN_SRC python
import numpy as np
ppar = [4, 3, -2, 10]
p = np.poly1d(ppar)
print p(3)
print np.polyval(ppar, 3)
x = 3
print 4*x**3 + 3*x**2 -2*x + 10
#+END_SRC
#+RESULTS:
: 139
: 139
: 139
numpy makes it easy to get the derivative and integral of a polynomial.
Consider: $y = 2x^2 - 1$. We know the derivative is $4x$. Here we compute the derivative and evaluate it at x=4.
#+BEGIN_SRC python
import numpy as np
p = np.poly1d([2, 0, -1])
p2 = np.polyder(p)
print p2
print p2(4)
#+END_SRC
#+RESULTS:
:
: 4 x
: 16
The integral of the previous polynomial is $\frac{2}{3} x^3 - x + c$. We assume $C=0$. Let us compute the integral $\int_2^4 2x^2 - 1 dx$.
#+BEGIN_SRC python
import numpy as np
p = np.poly1d([2, 0, -1])
p2 = np.polyint(p)
print p2
print p2(4) - p2(2)
#+END_SRC
#+RESULTS:
: 3
: 0.6667 x - 1 x
: 35.3333333333
One reason to use polynomials is the ease of finding all of the roots using numpy.roots.
#+BEGIN_SRC python
import numpy as np
print np.roots([2, 0, -1]) # roots are +- sqrt(2)
# note that imaginary roots exist, e.g. x^2 + 1 = 0 has two roots, +-i
p = np.poly1d([1, 0, 1])
print np.roots(p)
#+END_SRC
#+RESULTS:
: [ 0.70710678 -0.70710678]
: [ 0.+1.j 0.-1.j]
There are applications of polynomials in thermodynamics. The van der waal equation is a cubic polynomial $f(V) = V^3 - \frac{p n b + n R T}{p} V^2 + \frac{n^2 a}{p}V - \frac{n^3 a b}{p} = 0$, where $a$ and $b$ are constants, $p$ is the pressure, $R$ is the gas constant, $T$ is an absolute temperature and $n$ is the number of moles. The roots of this equation tell you the volume of the gas at those conditions.
#+BEGIN_SRC python
import numpy as np
# numerical values of the constants
a = 3.49e4
b = 1.45
p = 679.7 # pressure in psi
T = 683 # T in Rankine
n = 1.136 # lb-moles
R = 10.73 # ft^3 * psi /R / lb-mol
ppar = [1.0, -(p*n*b+n*R*T)/p, n**2*a/p, -n**3*a*b/p];
print np.roots(ppar)
#+END_SRC
#+RESULTS:
: [ 5.09432376+0.j 4.40066810+1.43502848j 4.40066810-1.43502848j]
Note that only one root is real (and even then, we have to interpet 0.j as not being imaginary. Also, in a cubic polynomial, there can only be two imaginary roots). In this case that means there is only one phase present.
*** Summary
Polynomials in numpy are even better than in Matlab, because you get a polynomial object that acts just like a function. Otherwise, they are functionally equivalent.