带数值积分的Sympy级数展开

Question

我想对函数F(e,Eo)进行一系列扩展，直到e的某个幂，然后在Eo变量上进行数值积分。我的想法是用SymPy在e中做幂级数，然后用MPMath进行Eo上的数值积分。

下面是一个示例代码。我收到的信息是，它不能从表达式中创建强积金。我猜问题与SymPy的系列在结尾有一个O(e**5)项有关，后来我希望数值积分显示一个e函数，而不是一个数字。

import sympy as sp
import numpy as np
from mpmath import *
e = sp.symbols('e')
Eo = sp.symbols('Eo')
expr = sp.sin(e-2*Eo).series(e, 0, 5)
F = lambda Eo : expr
I = quad(F, [0, 2*np.pi])
print(I)

很明显，我需要一种不同的方法，但我仍然需要对我的实际代码进行数值集成，因为它的表达式要复杂得多，无法进行解析集成。

编辑:我应该为示例代码选择一个由实变量组成的复杂函数，我正在尝试这样的函数(扩展和集成)，例如：

expr = (cos(Eo) - e - I*sqrt(1 - e ** 2)*sin(Eo)) ** 2 * (cos(2*(Eo - e*sin(Eo))) + I*sin(2*(Eo - e*sin(Eo))))/(1 - e*cos(Eo)) ** 4

Answer 1

下面是一种类似于Wrzlprmft的答案的方法，但通过SymPy的多项式模使用了一种不同的处理系数的方法：

from sympy import sin, pi, symbols, Integral, Poly

def integrate_coeff(coeff):
    partial_integral = coeff.integrate((Eo, 0, 2*pi))
    return partial_integral.n() if partial_integral.has(Integral) else partial_integral

e,Eo = symbols("e Eo")
expr = sin(e-sin(2*Eo))  
degree = 5

coeffs = Poly(expr.series(e, 0, degree).removeO(), e).all_coeffs()
new_coeffs = map(integrate_coeff, coeffs)
print((Poly(new_coeffs, e).as_expr() + e**degree).series(e, 0, degree))

主要代码有三行：(1)提取e到给定程度的系数；(2)在必要的情况下，对每个系数进行数值积分；(3)打印结果，将其表示为级数而不是多项式(因此增加e**度的技巧使SymPy知道该级数仍在继续)。输出：

-6.81273574401304e-108 + 4.80787886126883*e + 3.40636787200652e-108*e**2 - 0.801313143544804*e**3 - 2.12897992000408e-109*e**4 + O(e**5) 

Answer 2

我希望数值积分显示e的函数，而不是数字。

这在根本上是不可能的。您的集成要么是分析性的，要么是数字的，如果是数值的话，它只会为您处理和产生数字(数字和数字由于某种原因是相似的)。

如果您想将集成划分为数值和分析组件，那么您必须自己动手，- or希望SymPy能够根据需要自动拆分集成，不幸的是它还没有能力。我就是这样做的(详细信息在代码中有注释)：

from sympy import sin, pi, symbols, Integral
from itertools import islice

e,Eo = symbols("e Eo")
expr = sin(e-sin(2*Eo))

# Create a generator yielding the first five summands of the series.
# This avoids the O(e**5) term. 
series = islice(expr.series(e,0,None),5)

integral = 0
for power,summand in enumerate(series):
    # Remove the e from the expression
    Eo_part = summand/e**power
    # … and ensure that it worked:
    assert not Eo_part.has(e)

    # Integrate the Eo part:
    partial_integral = Eo_part.integrate((Eo,0,2*pi))

    # If the integral cannot be evaluated analytically, …
    if partial_integral.has(Integral):
        # … replace it by the numerical estimate:
        partial_integral = partial_integral.n()

    # Re-attach the e component and add it to the sum:
    integral += partial_integral*e**power

请注意，我根本没有使用NumPy或MPMath (尽管SymPy在数值估计中使用后者)。除非您真的知道自己在做什么，否则将这两者与SymPy混为一谈是个坏主意，因为他们甚至不知道SymPy表达式。