带有数值积分的 Sympy 级数展开

Question

我想对函数 F(e,Eo) 进行级数扩展，直到 e 的某个幂，并以数字方式对 Eo 变量进行积分。 我的想法是使用 SymPy 在e 中进行幂级数，然后使用 MPMath 对Eo 进行数值积分。

以下是示例代码。我收到无法从表达式创建 mpf 的消息。我想这个问题与 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

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

这基本上是不可能的。 您的积分要么是分析的，要么是数值的，如果它是数值的，它只会为您处理和生成数字（numerical 和 number 这两个词是相似的）。

如果您想将集成拆分为数值和分析组件，您必须自己进行 - 或者希望 SymPy 根据需要自动拆分集成，which it unfortunately is not yet capable of。 我就是这样做的（细节在代码中注释）：

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 表达式。

Answer 2

通过 SymPy 的polynomial module，这是一种类似于 Wrzlprmft 的答案但处理系数的方式不同的方法：

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**degree 的技巧，以使 SymPy 知道系列继续）。输出：

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