python：集成分段函数

Question

我想对一个已定义的函数进行分段积分，该函数乘以勒让德多项式。 不幸的是，我在documentation 中找不到如何使用 x 的第 n 个勒让德多项式。我想在n = 1,..., 50 时对x 的每个勒让德多项式进行积分，所以我设置了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

我想我需要定义另一个函数，比如说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)) 

所以我会返回一些 u(x)，前 50 个项通过勒让德多项式扩展我的分段函数。

编辑 1：

如果这不能完成，我可以使用罗德里格斯公式来计算第 n 个勒让德多项式。但是，当我在 Python 中寻找计算 n 次导数时，我找不到任何有用的东西。

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

所以如果有人知道如何在 Python 中实现这样的方案，这是一个选择。

编辑 2：

使用 Saullo Castro 的回答，我有：

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

但是，如果我打印 c，则值是错误的。值应为1.5, 0, -7/8, 0, ...

另外，当我绘制g时，我想做x = np.linspace(-1, 1, 500000)，所以情节很详细，但c只有50。如何实现？

Answer 1

如果我正确理解您的问题，您想计算 f(x) * Ln(x) 的积分，其中 f(x) 是您使用 python 函数定义的分段函数。我假设您对这个特定的步进函数并不特别感兴趣。

您可以使用 legval 和系数参数的单位矩阵来获取勒让德多项式的值。

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)

然后您可以执行积分积分。使用 gauss-legendre 求积可能更有效，因为 Legendre 多项式的积分对于 Ln(x) 将是精确的，其中 n 小于求积大小。

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, ... 我认为这是您正在寻找的结果。

Answer 2

您可以像这样进行集成：

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

这给出了：

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