python: integrating a piecewise function

Question

I want to integrate a piecewise a defined function that is multiplied by the Legendre polynomials. Unfortunately, I can't find how to use the nth Legendre polynomial of x in the documentation. I want to integrate each Legendre polynomial of x when n = 1,..., 50 so I have set n = np.arange(1, 51, 1).

import numpy as np
import pylab
from scipy import integrate

n = np.arange(1, 51, 1)                                                   


def f(x):
    if 0 <= x <= 1:
        return 1
    if -1 <= x <= 0:
        return -1

I suppose I need to define another function let's say u(x).

c = []


def u(x):
    c.append((2. * n + 1) / 2. * integrate.quad(f(x) * insert Legendre polynomials here, -1., 1.)) 
    return sum(c * Legendre poly, for nn in range(1, 51)) 

So I would then return some u(x) with the first 50 terms expanding my piecewise function by Legendre polynomials.

Edit 1:

If this can't be done, I could use Rodrigues's Formula to compute the nth Legendre polynomial. However, I couldn't find anything useful when I was looking for computing nth derivatives in Python.

P_n(x) = \frac{1}{2^n n!}\frac{d^n}{dx^n}(x^2 - 1)^n

So this is an option if someone knows how to implement such a scheme in Python.

Edit 2:

Using Saullo Castro's answer, I have:

import numpy as np
from scipy.integrate import quad

def f(x, coef):
    global p
    p = np.polynomial.legendre.Legendre(coef=coef)
    if 0 <= x <= 1:
        return 1*p(x)
    if -1 <= x <= 0:
        return -1*p(x)

c = []
for n in range(1, 51):
    c.append((2. * n + 1.) / 2. * quad(f, -1, 1, args=range(1,n+1))[0])

def g(x)
    return sum(c * p(x) for n in range(1, 51))

However, if I print c, the values are wrong. The values should be 1.5, 0, -7/8, 0, ...

Also, when I plot g, I would like to do x = np.linspace(-1, 1, 500000) so the plot is detailed but c is only 50. How can this be achieved?

Answer 1

If I understand your question correctly, you want to calculate the integral of f(x) * Ln(x) where f(x) is a piecewise function you're defining with a python function. I'm assuming you're not specifically interested in this particular step function.

You can get the values of the Legendre polynomials using legval with the identity matrix for the coefficient argument.

import numpy as np
import matplotlib

x = np.linspace(-1, 1, 201)

L = np.polynomial.legendre.legval(x, np.identity(50))

plt.plot(x, L.T)

You can then perform the integral with quadrature. Using gauss-legendre quadrature might be more efficient since the integral of a legendre polynomial will be exact for Ln(x) where n is less than the quadrature size.

import numpy as np    
from numpy.polynomial.legendre import leggauss, legval

def f(x):
    if 0 <= x <= 1:
        return 1
    if -1 <= x <= 0:
        return -1

# of course you could write a vectorized version of
# this particular f(x), but I assume you have a more
# general piecewise function
f = np.vectorize(f)

deg = 100
x, w = leggauss(deg) # len(x) == 100

L = np.polynomial.legendre.legval(x, np.identity(deg))
# Sum L(xi)*f(xi)*wi
integral = (L*(f(x)*w)[None,:]).sum(axis=1)

c = (np.arange(1,51) + 0.5) * integral[1:51]

x_fine = np.linspace(-1, 1, 2001) # 2001 points
Lfine = np.polynomial.legendre.legval(x_fine, np.identity(51))

# sum_1_50 of c(n) * Ln(x_fine)
cLn_sum = (c[:,None] * Lfine[1:51,:]).sum(axis=0)

c = 1.5, 0, -8.75e-1, 0, ... which I think is the result you're looking for.

Answer 2

You can do your integration like:

import numpy as np
from scipy.integrate import quad

def f(x, coef):
    n = coef[-1]
    p = (2*n+1)/2.*np.polynomial.legendre.Legendre(coef=coef)
    if 0 <= x <= 1:
        return 1*p(x)
    if -1 <= x <= 0:
        return -1*p(x)

c = []
for n in range(1, 51):
    c.append(quad(f, -1, 1, args=range(1,n+1))[0])

And this gives:

print c
#[0.0, 5.0, 6.9999999999999991, 4.5, ... , 62.975635570615466, 64.274102283412574, 77.143785770271251]