Symbolic integration slow with SymPy
Asked Answered
O

1

6

I'm not a SymPy expert, but I've used it successfully in some of my lectures in the last years. However, sometimes it seems to be very slow with symbolic integration. Here's an example that Mathematica computes almost instantly while SymPy needs a long time (more than half a minute) for it on my machine.

from sympy import *
x = symbols("x")

def coeff (f, k, var = x):
    return integrate(f(var) * cos(k * var), (var, -pi, pi)) / pi

f = lambda x: (11*sin(x) + 6*sin(2*x) + 2*sin(3*x))/10

[coeff(f, k) for k in range(0, 5)]

Am I doing something wrong or is this the expected behavior? Is there some trick to speed this up?

SymPy version is 1.0, Python is 3.5.1 64-Bit (Anaconda) on Windows.

Orling answered 10/12, 2017 at 14:0 Comment(2)
What does slow mean in this case?Sturm
What I wrote in my question. More than half a minute in SymPy versus instant reply on Mathematica. Is that not enough information?Orling
P
5

You can help guide the integration by telling SymPy to expand products of sin-cos to sums (see also [1],[2]) before doing the integration:

For example,

In [58]: fprod(f, 4)
Out[58]: (11*sin(x)/10 + 3*sin(2*x)/5 + sin(3*x)/5)*cos(4*x)

In [59]: FU['TR8'](fprod(f, 4))
Out[59]: -sin(x)/10 - 3*sin(2*x)/10 - 11*sin(3*x)/20 + 11*sin(5*x)/20 + 3*sin(6*x)/10 + sin(7*x)/10

The integration is simpler in this form.


Therefore you could use:

import sympy as sym
x = sym.symbols("x")


def f(x):
    return (11*sym.sin(x) + 6*sym.sin(2*x) + 2*sym.sin(3*x))/10

def fprod(f, k, var=x):
    return f(var) * sym.cos(k * var)

FU = sym.FU
def coeff (f, k, var=x):
    return sym.integrate(FU['TR8'](fprod(f, k)), (var, -sym.pi, sym.pi)) / sym.pi

[coeff(f, k) for k in range(0, 5)]

Here is a benchmark using FU['TR8']:

In [52]: %timeit [coeff(f, k) for k in range(0, 5)]
10 loops, best of 3: 78.8 ms per loop

Using the original code (without FU['TR8']):

In [54]: %timeit [coeff(f, k) for k in range(0, 5)]
1 loop, best of 3: 19.8 s per loop
Palindrome answered 10/12, 2017 at 14:33 Comment(1)
Thanks, that's what I was looking for.Orling

© 2022 - 2024 — McMap. All rights reserved.