Python 等效于 MATLAB 的 Legendre 函数答案

【问题标题】：Python equivalent of MATLAB's Legendre functionPython 等效于 MATLAB 的 Legendre 函数
【发布时间】：2014-11-21 12:16:49
【问题描述】：

目前，我正在尝试使用 Python 分析时间序列数据。作为这样做的指导方针，我将自己定位在一个 MATLAB 脚本上，它几乎可以完成我想做的所有事情。到目前为止它运行良好，但现在我遇到了该脚本中使用的勒让德多项式。

我尝试了其中的NumPy implementation，但我找不到（或多或少）产生与the MATLAB function 相同结果的方法。

基本上，这就是我想知道的。如何让我的 Python 代码给出与 MATLAB 代码相同的结果？

作为一个小示范，

    k= [0 1 1;1 1 0 ;0 0 1]
    legendre(2,k)

给予：

ans(:,:,1) =

-0.5000    1.0000   -0.5000
0         0         0
3.0000         0    3.0000


ans(:,:,2) =

1.0000    1.0000   -0.5000
     0         0         0
     0         0    3.0000


ans(:,:,3) =

1.0000   -0.5000    1.0000
     0         0         0
     0    3.0000         0

而我的 Python 版本是这样的：我尝试的方式是这样的：

    legendre = np.polynomial.legendre.Legendre([0,1,2])
    legendre(k)

而产量：

   array([[-1.,  3.,  3.],
   [ 3.,  3., -1.],
   [-1., -1.,  3.]])

我看到一些有点奇怪的东西，但不幸的是我不知道如何测试它们，因为这是我第一次听说像勒让德多项式这样的东西，而且 NumPy 的文档和维基百科都不是对理解它有很大帮助。

【问题讨论】：

标签： python matlab numpy polynomials

【解决方案1】：

import numpy as np
from scipy.special import lpmv

def legendre(deg,x):
    return np.asarray([lpmv(i,deg,x) for i in range(deg+1)])

x=np.array([[0,1,1],[1,1,0],[0,0,1]])
legendre(2,x)

它给了你想要的。

【讨论】：

【解决方案2】：

我也遇到了这个问题。以这个问题为起点，提出了以下问题。请注意：我正在使用 MATLAB 函数，如下所示：

legendre(10,linspace(-1,1,10))

我需要在 Python 中生成等价物。代码如下：

import numpy as np
from scipy import special as sp

def legendre(N,X) :
    matrixReturn = np.zeros((N+1,X.shape[0]))
    for i in enumerate(X) :
        currValues = sp.lpmn(N,N,i[1])
        matrixReturn[:,i[0]] = np.array([j[N] for j in currValues[0]])
    return matrixReturn

我对 Python 很陌生，所以我确信上述内容可以改进。

【讨论】：

【解决方案3】：

好的，我认为使用此模块复制这些结果将有困难，因为从名称判断仅处理勒让德多项式（这些是勒让德方程的解，其中 mu = 0，也称为 0 阶解)

我不知道matlab，但是看文档，你的输入是计算legendre函数的结果，最多指定度数。

在 python 中，您似乎正在做的是创建第零一阶和二阶 Legendre 多项式的组合

0*l_0 + 1*l_1 + 2*l_2

您可以在指定的点计算 Legendre 多项式：

l0 = np.polynomial.legendre.Legendre([0,1])

你可以验证一下

l0(0.5) == 0.5

我希望这很有用 - 请随时提出更多问题

编辑：

def coefficients(order):
    for i in range(1, order):
         base = np.zeros(order)
         base[i] = 1
         yield base

def a(n, m):
    return 1.0*(2*n+1) / ((n*(n+1))**m)

def g(const_dist, m, order):
     legendres = [np.polynomial.legendre.Legendre(n) for n in coefficients(order)]
     terms = [a(n+1,m)*legendres[n](const_dist) for n,_ in enumerate(legendres)]
     return sum(terms)


>>> g(0.4, 4, 6)
0.073845698737654328

我希望这对你有用，如果我搞砸了，请告诉我

【讨论】：

谢谢，我现在明白了一点。那么，为了在python中模仿matlab的legendre，我自己不能使用这个模块，但必须编写自己的legendre函数？您认为将 Rodrigues 公式（参见 wiki）放入 python 中是否足够（并且可行），以便我可以分别创建多项式并在指定点评估它们中的每一个？顺便提一句。我实际上必须翻译的公式是pdf 幻灯片 14 上的第一个公式。任何想法如何做到这一点？
实际上，我认为 python 模块在这里就足够了，因为您似乎只是在计算legendre 多项式，而不是legendre 函数。我将尝试将此函数的可能计算添加到我的答案中
您选择“order”为“6”而“m”为“4”有什么原因吗？正常情况是不是你把每一个多项式都取到最后一个，或者我什么也得不到？
忘记我的最后评论，这是一个愚蠢的问题。我刚看了一下实际的公式，发现“m”是完全不相关的。对不起

【解决方案4】：

@user3684792 感谢您的代码，但这并不是我们真正需要的，例如cosdist 通常是一个矩阵，所以 sum(terms) 是不够的（虽然很容易修复）。

根据您的评论和this 罗格朗德多项式的定义，我自己尝试了。我最终得到的是这段代码。请问您对此有何看法？

    def P(n,x):
        if n == 0:
            return  1
        elif n==1:
            return  x
        elif n>1:
            return  (2*n-1)/n * x * P(n-1,x) - (n-1)/n * P(n-2,x)

    #some example data
    order = 4
    cosdist= np.array(((0.4,0.1),(-0.2,0.3)))
    m = 3
    dim1_cosdist, dim2_cosdist = cosdist.shape

    Gf = np.zeros((order, dim1_cosdist, dim2_cosdist))
    for n in range(1,order):
        Gf[n] = 1.0*(2*n+1) / ((n*(n+1))**m) * P(n,cosdist) 

    G = np.sum(Gf,axis = 0)

如果 cosdist 只是一个整数，此脚本会给出与您相同的结果。令我恼火的是，这些结果仍然与 Matlab 代码产生的结果有些不同，即结果数组甚至具有不同的维度。谢谢。编辑：意外地，我将m 与order 混淆了。现在应该是正确的

【讨论】：

很好地使用了递归关系。这是获取它们的超级简单方法，可能更适合您使用。（虽然它可能会更慢，所以如果速度很重要，您应该对此进行分析）。
我同意用 codistance 作为矩阵来计算它可能会更好，但假设矩阵只是一组系数（即数字），您实际上可以通过运行方程来计算 Gij每个距离 Ij 数。
wrt 维度，您可以通过确保 Gij 具有与 Codistance ij 相同的维度来进行完整性检查。您注意到尺寸与您的 matlab 代码不同。我记得当我第一次查看该函数时，看到它还计算 Legendre 函数，所以这可能是您在那里也看到的（如果是这样，那不是您想看到的）。
计算时间并不重要。我必须处理的数组大小为 128 x 128，并且只计算一次，但我会记住这一点，谢谢。顺便提一句。我应该以某种方式结束这个话题吗？或者现在处理它的最佳方法是什么，因为我已经在这里完成并且没有回答实际（原始）问题？
我想如果您认为原始 q 没有得到回答，请保持打开状态，也许它会对其他人有所帮助。这个问题的答案是，我认为 matlab 也提供了 legendre fn 解决方案，这不是你想要的

【解决方案5】：

SciPy 具有associated Legendre polynomials。它与 MATLAB 版本不同，但它应该提供您想要的大部分功能。

【讨论】：

【解决方案6】：

我遇到了同样的问题，并成功构建了以下内容：

from scipy import special

def legendre(n,X) :
res = []
for m in range(n+1):
    res.append(special.lpmv(m,n,X))
return res

对于我的应用程序来说，这非常有效 - 也许你们中的一些人也可以使用它。

【讨论】：