Sympy系列扩展与数值整合

I want to make a series expansion for a function F(e,Eo) up to a certain power of e and integrate over the Eo variable numerically. What I thought was using SymPy to make the power series in e, and then use MPMath for the numerical integration over Eo.

我想对函数F（e，Eo）进行一系列扩展，直到e的某个幂，并在数字上整合Eo变量。我认为使用SymPy在e中制作幂级数，然后使用MPMath进行Eo上的数值积分。

Below is an example code. I receive the message that it can not create mpf from the expression. I guess the problem has to do with the fact that with the series from SymPy has an O(e**5) term at the end, and later that I want the numerical integration to show a function of e instead of a number.

下面是一个示例代码。我收到的消息是它无法从表达式创建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)

It’s evident that I need a different approach, but I would still need the numerical integration for my actual code because it has much more complicated expressions that could not be integrated analytically.

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

Edit: I should have chosen a complex function of real variables for the example code, I am trying this (the expansion and integration) for functions such as:

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

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

2 个解决方案

#1

Here is an approach similar to Wrzlprmft's answer but with a different way of handling coefficients, via SymPy's polynomial module:

这是一种类似于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))

The main code is three lines: (1) extract coefficients of e up to given degree; (2) integrate each, numerically if we must; (3) print the result, presenting it as a series rather than a polynomial (hence the trick with adding e**degree, to make SymPy aware that the series continues). Output:

主要代码是三行：（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)

#2

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

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

This is fundamentally impossible. Either your integration is analytical or numerical, and if it is numerical it will only handle and yield numbers for you (the words numerical and number are similar for a reason).

这基本上是不可能的。你的集成是分析的还是数字的，如果它是数字的，它只会为你处理和产生数字（单词数字和数字是相似的）。

If you want to split the integration into numerical and analytical components, you have to do so yourself – or hope that SymPy automatically splits the integration as needed, which it unfortunately is not yet capable of. This is how I would do it (details are commented in the code):

如果你想将集成拆分为数值和分析组件，你必须自己这样做 - 或者希望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

Note that I did not use NumPy or MPMath at all (though SymPy uses the latter under the hood for numerical estimates). Unless you really really know what you’re doing, mixing those two with SymPy is a bad idea as they are not even aware of SymPy expressions.

请注意，我根本没有使用NumPy或MPMath（尽管SymPy使用后者进行数值估算）。除非你真的知道你在做什么，否则将这两者与SymPy混合是一个坏主意，因为他们甚至不知道SymPy表达式。

#1