0

I'm trying to use Python to calculate the Rodrigues formula, P_n(x).

http://en.wikipedia.org/wiki/Rodrigues%27_formula

That is, I would like a function which takes into two input parameters, n and x, and returns the output of this formula.

However, I don't think SciPy has this function yet. SpiPy does offer a Legendre module:

http://docs.scipy.org/doc/numpy/reference/routines.polynomials.legendre.html

I don't think any of these is the Rodrigues formula. Am I wrong?

Is there a standard way SciPy offers to do this?

EDIT: I would like the input parameters to be arrays, not just single input values.

1

1 Answer 1

4

If you simply want P_n(x), then you can create a suitable object representing the P_n polynomial using scipy.special.legendre and call it with your values of x:

In [1]: from scipy.special import legendre
In [2]: n = 3
In [3]: Pn = legendre(n)
In [4]: Pn(2.5)
Out[4]: 35.3125        # P_3(2.5)

The object Pn is, in a sense, the "output" of the Rodrigues formula: it is a polynomial of the required order, which can be evaluated at a provided value of x. If you want a single function that takes n and x, you can use eval_legendre:

In [5]: from scipy.special import eval_legendre
In [6]: eval_legendre(3, 2.5)
Out[6]: 35.3125

As noted in the docs, this is the recommended way to do it for large-ish n (e.g. n > 20), instead of creating a polynomial object with all the coefficients which does not handle rounding errors and numerical stability as well.

EDIT: Both approaches work with arrays (at least for the x argument). For example:

In [7]: x = np.array([0, 1, 2, 5, 10])
In [8]: Pn(x)
Out[8]: 
array([  0.00000000e+00,   1.00000000e+00,   1.70000000e+01,
     3.05000000e+02,   2.48500000e+03])
5
  • It would be nice if I could get something that works with n as an array as well. This is a great start though. Thanks! Apr 20, 2015 at 20:12
  • scipy.special.eval_legendre
    – pv.
    Apr 20, 2015 at 21:45
  • Thanks, @pv: I'd forgotten about the eval_ methods: edited.
    – xnx
    Apr 21, 2015 at 0:51
  • If you only need function values, and not the power basis representation and roots etc., using eval_legendre is recommended (try legendre(1000)(0.1) and eval_legendre(1000, 0.1)). Then read the warning notice in the documentation, and look at x=0.1*legendre(50); print(x(0.9), 0.1*legendre(50)(0.9), 0.1*eval_legendre(50, 0.9))
    – pv.
    Apr 21, 2015 at 8:20
  • This warning is well-made, though I don't know whether it applies to OP's use-case. I've edited in a few comments, but if you have anything to add feel free to edit or add an answer yourself.
    – xnx
    Apr 21, 2015 at 15:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.