如何用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算法，在某些情况下它的速度非常慢。还有一种“手动”积分方法，它不像Risch那样完整，但偶尔会出现极端缓慢的情况。这里有一些描述： https://docs.sympy.org/latest/modules/integrals/integrals.html#internals

您给出的示例中的问题是，它被困在heurisch中。所以让我们试试“手动”吧：

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留下了一个未赋值的积分。我们可以使用普通的integrate通过调用doit来完成：

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

虽然这个答案看起来与前一个不同，但是有无数种使用trig恒等式重写trig表达式的方法，但是它们应该等价于一个加法常数（对于反导数是必需的）。你知道吗