CODEWITHSUNDEEP
pythonmatlabnumpypolynomials
userE
userE

Reputation: 249

Python equivalent of MATLAB's Legendre function

Currently, I am trying to analyze time series data with Python. As a guideline for doing so, I oriented myself on a MATLAB script that does pretty much everything I'd like to do. So far it worked fine, but now I bumped into this Legendre polynomial that was used in that script.

I tried the NumPy implementation of it, but I couldn't find a way that (more or less) yielded the same results as the MATLAB function.

Basically, this is what I'd like to know. How can I make my Python code give the same results as the MATLAB code?

As a small demonstration,

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

gives:

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

Whereas my Python version of it goes like this: the way I tried it goes like so:

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

And yields:

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

I see a few things that are a bit weird, but unfortunately I have no clue about how to test for them, because this is the first time I heard of such a thing like Legendre polynomial and neither NumPy's documentation nor Wikipedia are a big help understanding it.

Upvotes: 2

Views: 4929

Answers (6)

flumer
flumer

Reputation: 176

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)

it gives what you want.

Upvotes: 2

maxm59
maxm59

Reputation: 46

I had the same problem and was successful with building the following:

from scipy import special

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

For my applications this works very well - maybe some of you can use this too.

Upvotes: 3

J Pierret
J Pierret

Reputation: 61

I just ran into this problem as well. Used this question as a starting point and came up with the following. Please note: I'm using the MATLAB function like so:

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

And I needed to generate the equivalent in Python. Here's the code:

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

I'm very new to Python, so I'm sure the above could be improved.

Upvotes: 4

Charles Harris
Charles Harris

Reputation: 984

SciPy has the associated Legendre polynomials. It isn't the same as the MATLAB version, but it should provide most of what you want.

Upvotes: 2

userE
userE

Reputation: 249

@user3684792 Thanks for the code, but this is not exactly what is needed, e.g. cosdist is normally a matrix, so that sum(terms) would not suffice (easy fix though).

Based on your comment and this definition of Legrande polynomials, I tried it myself. What I ended up with is this code. May I ask for your opinion about it?

    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)

In case of cosdist is just an integer, this script gives the same results as yours. What irritates me, is that these results are still somewhat different than those resulting from the Matlab code, i.e. the resulting array has even different dimensions. Thanks. Edit: By accident, I confused m with order. Now it should be correct

Upvotes: 2

user3684792
user3684792

Reputation: 2611

Ok, well I think you will have trouble replicating these reslults using this module, as judging by the name only deals with the legendre polynomials (These are the solutions to the legendre equation where mu = 0, otherwise known as order 0 solutions)

I don't know matlab, but looking at the documentation, your input is calculating the results of the legendre functions of up to the order of the degree specified.

In python, what you seem to be doing is creating a composition of the zeroeth first and second order legendre polynomials

0*l_0 + 1*l_1 + 2*l_2

you can evaluate the legendre polynomials at the points specified:

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

and you can verify that 

l0(0.5) == 0.5

I hope this is useful - feel free to ask more questions

Edit:

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

I hope this works for you, let me know if I have messed something up

Upvotes: 2

Related Questions