python如何对大三角方程进行积分？

Question

在我的程序中，我必须重新创建 galerkin-method，我必须用它做一些我现在用 sympy.integrate 做的集成。

我的问题是，在方程包含多个三角函数或添加的 e 函数（它应该能够解决所有问题）的情况下，sympy 会永远计算。我也在使用 sympy.symbols，因为通过这些集成，我正在创建一个必须求解的方程组。我需要一个反导数，而不仅仅是一个值的解决方案。

是否有数值积分方法或其他返回准确值的方法？我尝试使用 sympy.Integral(equation).evalf(1)，但错误太高或返回的十进制数太长，这让我回到了运行时间太长的状态。

函数可以如下所示：

(-4*x**3*sin(2*x) + 2*cos(x)*cos(2*x))*(cos(x) + 1)

我必须整合其中的 20 个。

Answer 1

如果您不需要符号积分，并且您知道积分界限，那么您可以进行数值积分。一个务实的起点是梯形法则。 https://en.wikipedia.org/wiki/Trapezoidal_rule。可以通过更精细的步骤（在范围内）任意提高准确性。

通过使用更高的阶数来提高准确性在编程上更复杂一些，但在数值上更有效。一个实现通常只有在计算时间高于编程时间时才开始有意义）

import math
import matplotlib

sin=math.sin  
cos=math.cos

def f(x):
    return (-4*x**3*sin(2*x) + 2*cos(x)*cos(2*x))*(cos(x) + 1)

x0 = 0              # lower bound
x1 = 2.8            # upper bound (not exactly resolved in this example)
dx = 0.01           # integration step size

F_int = 0
x=x0

while x < x1:
    f_n = f(x)
    f_n_plusdx = f(x+dx)
    dF = 0.5*(f_n+f_n_plusdx)*dx     # calc the trapezoid area
    F_int += dF                      # sum of the integration
    x+=dx                            # next dx.

print(F_int)

Answer 2

SymPy 的集成函数尝试了多种不同的集成方法，其中一种是 Risch 算法，它在某些情况下可能非常慢。还有 "manual" 积分方法，它不如 Risch 完整，但偶尔出现极度缓慢的情况较少。这里有一些描述： https://docs.sympy.org/latest/modules/integrals/integrals.html#internals

您给出的示例中的问题是它卡在启发式中。所以让我们尝试 "manual" 代替：

In [1]: expr = (-4*x**3*sin(2*x) + 2*cos(x)*cos(2*x))*(cos(x) + 1)                                                                

In [2]: expr                                                                                                                      
Out[2]: 
⎛     3                             ⎞             
⎝- 4⋅x ⋅sin(2⋅x) + 2⋅cos(x)⋅cos(2⋅x)⎠⋅(cos(x) + 1)

In [3]: %time anti1 = expr.integrate(x, manual=True)                                                                              
CPU times: user 39.7 s, sys: 232 ms, total: 39.9 s
Wall time: 43.1 s

In [4]: anti1                                                                                                                     
Out[4]: 
   3    3                                                              3                                           ⌠            
8⋅x ⋅cos (x)      3               2                           x   4⋅sin (x)                           sin(4⋅x)     ⎮  2    3    
──────────── + 2⋅x ⋅cos(2⋅x) - 3⋅x ⋅sin(2⋅x) - 3⋅x⋅cos(2⋅x) + ─ - ───────── + 2⋅sin(x) + 2⋅sin(2⋅x) + ──────── - 8⋅⎮ x ⋅cos (x) dx
     3                                                        2       3                                  8         ⌡

所以这花了 40 秒，但结果没有完全积分：manualintegrate 在那里留下了一个未计算的积分。我们可以通过调用 doit:

使用正常的 integrate 来完成它

In [5]: %time anti1.doit()                                                                                                        
CPU times: user 4.46 s, sys: 142 ms, total: 4.61 s
Wall time: 4.81 s
Out[5]: 
   3    3                          2    3                                                    2                      3           
8⋅x ⋅cos (x)      3            16⋅x ⋅sin (x)      2           2         2            32⋅x⋅sin (x)⋅cos(x)   112⋅x⋅cos (x)        
──────────── + 2⋅x ⋅cos(2⋅x) - ───────────── - 8⋅x ⋅sin(x)⋅cos (x) - 3⋅x ⋅sin(2⋅x) - ─────────────────── - ───────────── - 3⋅x⋅c
     3                               3                                                        3                  9              

                     3                    2                                      
          x   284⋅sin (x)   112⋅sin(x)⋅cos (x)                           sin(4⋅x)
os(2⋅x) + ─ + ─────────── + ────────────────── + 2⋅sin(x) + 2⋅sin(2⋅x) + ────────
          2        27               9                                       8

所以又花了几秒钟才得到那个结果。这现在是一个完整的反导数，我们可以验证：

In [6]: simplify(expr - _.diff(x))                                                                                                
Out[6]: 0

这意味着我们可以用 expr.integrate(x, manual=True).doit() 在大约 50 秒内完成这个特定的积分。

实际上，如果将此特定示例从 sin/cos 重写为 exp:

，则可以在 5 秒内完成

In [1]: expr = (-4*x**3*sin(2*x) + 2*cos(x)*cos(2*x))*(cos(x) + 1)                                                                

In [2]: %time expr.rewrite(exp).expand().integrate(x).expand().rewrite(sin).simplify()                                            
CPU times: user 5.3 s, sys: 21.2 ms, total: 5.32 s
Wall time: 5.33 s
Out[2]: 
                                 3                                             2                                                
   3             3            2⋅x ⋅cos(3⋅x)      2             2            2⋅x ⋅sin(3⋅x)                                4⋅x⋅cos
2⋅x ⋅cos(x) + 2⋅x ⋅cos(2⋅x) + ───────────── - 6⋅x ⋅sin(x) - 3⋅x ⋅sin(2⋅x) - ───────────── - 12⋅x⋅cos(x) - 3⋅x⋅cos(2⋅x) - ───────
                                    3                                             3                                           9 


(3⋅x)   x                            13⋅sin(3⋅x)   sin(4⋅x)
───── + ─ + 13⋅sin(x) + 2⋅sin(2⋅x) + ─────────── + ────────
        2                                 27          8    

In [3]: simplify(expr - _.diff(x))                                                                                                
Out[3]: 0

尽管这个答案看起来与前一个答案不同，但使用三角恒等式重写三角表达式的方法有无数种，但它们应该等价于加法常数（根据反导数的要求）。

python如何对大三角方程进行积分？

How to integrate large trigonometric equations with python?

python

integration

sympy