使用Python SciPy来计算Rodrigues公式P_n(x)(勒让德多项式)

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

我正在尝试使用Python来计算Rodrigues公式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.

也就是说，我想要一个包含两个输入参数n和x的函数，并返回这个公式的输出。

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

然而，我认为SciPy还没有这个功能。SpiPy提供了一个勒让德模块:

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?

SciPy提供了一种标准的方法吗?

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

编辑:我希望输入参数是数组，而不是单个的输入值。

1 个解决方案

#1

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:

如果您只想要P_n(x)，那么可以使用scipy.special创建一个表示P_n多项式的合适对象。勒让德，用你的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:

从某种意义上说，对象Pn是Rodrigues公式的“输出”:它是所需顺序的多项式，可以在给定的x值上求值。

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.

正如在文档中所指出的，这是为大的n(例如n > 20)所推荐的方法，而不是创建一个具有所有系数的多项式对象，这些系数不处理舍入误差和数值稳定性。

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

编辑:这两种方法都可以使用数组(至少对于x参数)。例如:

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])

#1