具有数值积分的 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

I want the numerical integration to show a function of e instead of a number.

这根本不可能。您的积分要么是分析积分，要么是数值积分，如果是数值积分，它只会为您处理和生成数字（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

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

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)

具有数值积分的 Sympy 级数扩展

Sympy series expansion with numerical integration

python

sympy

numerical-integration