Sympy級數展開與數值積分

我想爲一個功能F(e,Eo)一系列拓展最高達e一定的功率和集成在Eo變量數值。我的想法是用SymPy使在e電源系列，然後在Eo使用MPMath的數值積分。Sympy級數展開與數值積分

下面是一個例子代碼。我收到它無法從表達式創建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

來源

2017-09-30 transistorNPN

旁註1：行'F =拉姆達EO：你想要什麼expr'幾乎肯定不會做，因爲'lambda表達式Eo'不一樣'expr中使用的一個'。你可能想用'F = lambda x：expr.subs（Eo，x）'或SymPy的'lambdify'。 – Wrzlprmft

旁註2：爲什麼使用MPMath的各種各樣的問題的人呢？它適用於任意精度的數學運算，這一切都是正確的，但也是大多數問題不需要的。 – Wrzlprmft

這裏是類似於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）提取到給定程度Ë向上的係數; （2）如果我們必須把每一個都進行數字化; （3）打印結果，其顯示爲一個系列，而不是一個多項式（因此招用新增電子郵件**度，使SymPy知道該系列繼續）。輸出：

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

來源

2017-09-30 17:42:24 FTP

不錯的解決方案，但是如何優化它，以免一旦該值看起來很小就不會浪費時間計算係數？就像您在代碼中使用的示例函數一樣，它會將係數與e-100訂單的結果進行整合。 – transistorNPN

我想要的數值積分顯示e而不是數字的功能。

這是根本不可能的。要麼你的積分是分析或數字，如果它是數字，它只會爲你處理和產生數字（數字和數字是類似的原因）。

如果您想要將集成拆分爲數字和分析組件，您必須自己做 - 或者希望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混合是一個壞主意，因爲他們甚至不知道SymPy表達式。

來源

2017-09-30 08:29:24 Wrzlprmft

我想知道它使用什麼方法對.integrate（）。n（）的數值估計。 – transistorNPN

Sympy級數展開與數值積分

回答

相關問題