python:集成分段函数 [英] python: integrating a piecewise function

问题描述

我想对乘以勒让德多项式的定义函数进行分段积分.不幸的是，我在 文档.我想在 n = 1,..., 50 时对 x 的每个勒让德多项式进行积分，所以我设置了 n = np.arange(1, 51, 1).

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

我想我需要定义另一个函数，比如 u(x).

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

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

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

编辑 1:

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

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

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

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

编辑 2:

使用 Saullo Castro 的回答，我有:

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

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

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

另外，当我绘制 g 时，我想做 x = np.linspace(-1, 1, 500000) 所以情节很详细但是 c 才50，怎么实现?

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?

推荐答案

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

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.

您可以使用 legval 与系数参数的单位矩阵.

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)

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

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

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

这篇关于python:集成分段函数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！