Python: Integrating A Piecewise Function

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 do

Solution 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
importmatplotlibx= 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

deff(x):
    if0 <= x <= 1:
        return1if -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.

Solution 2:

You can do your integration like:

import numpy as np
from scipy.integrate import quad

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

c = []
for n inrange(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]