How to optimise the numerical evaluation of a SymPy integral?

user_4955732 0
Viewed 178

I'm somewhat of a newbie to SymPy and was hoping someone could point out ways to optimise my code.

I need to numerically evaluate a somewhat involved expression with very high decimal places (150–300), and it is taking 30 seconds or longer per parameter set – which is very long given the parameter space to be calculated.

I have used lambdify with the mpmath backend and meijerg=True in the integral handling and it brought down run-times significantly. Are there any other methods that could be used? Ideally it would be great to push evaluation times below 1 second. My code is:

import mpmath
from mpmath import mpf, mp
mp.dps = 150 # ideally would like to have this set to 300 
import numpy as np
from sympy import besselj, symbols, hankel2, legendre, sin, cos, tan, summation, I
from sympy import  lambdify, expand, Integral
import time
x, alpha, k, m,n, r1, R, theta = symbols('x alpha k m n r1 R theta')

r1 = (R*cos(alpha))/cos(theta) #

Imn_part1 = (n*hankel2(n-1,k*r1)-(n+1)*hankel2(n+1,k*r1))*legendre(n, cos(theta))*cos(theta)
Imn_part2 = n*(n+1)*hankel2(n, k*r1)*(legendre(n-1, cos(theta)-legendre(n+1, cos(theta))))/k*r1
Imn_parts = expand(Imn_part1+Imn_part2)
Imn_expr = expand(Imn_parts*legendre(m,cos(theta))*(r1**2/R**2)*tan(theta))
Imn = Integral(Imn_expr, (theta, 0, alpha)).doit(meijerg=True)

# the lambdified expression
Imn_lambdify = lambdify([m,n,k,R,alpha], Imn,'mpmath')

When giving numerical inputs to the function – it takes a long time (30 s – 40 s).

substitute_dict = {'alpha':mpf(np.radians(10)), 'k':5,'R':mpf(0.1), 'm':20,'n':10}

print('starting calculation...')
start = time.time()
output = Imn_lambdify(substitute_dict['m'],
             substitute_dict['n'],
             substitute_dict['k'],
             substitute_dict['R'],
            substitute_dict['alpha'])
print(time.time()-start)

OS/package versions used:

  • Linux Mint 19.2
  • Python 3.8.5
  • SymPy 1.7.1
  • MPMath 1.2.1
mpmath sympy symbolic-math optimization python
1 Answers

Setting meijerg=True has just caused SymPy to not try as hard in evaluating the integral. It still can't evaluate it, but it has split it into 5 sub-integrals, which you can see if you print Imn. You might as well just leave it as one integral (leave off the doit()):

Imn = Integral(Imn_expr, (theta, 0, alpha))

For me, the split integral evaluates a little faster, but this is also about the same speed

Imn = Integral(simplify(Imn_expr), (theta, 0, alpha))

Ultimately, the thing that makes things slow is the number of digits that you are using. If you don't actually need these many digits, you shouldn't use them. Note that mpmath will automatically increase the precision internally to avoid cancellation, so it is unnecessary to do so yourself. I get the same value (with fewer digits) with the default dps of 15 as 150.

You can try substituting your values directly into your expression, if they do not change, and seeing if SymPy can simplify Imn_expr further with them.

As an aside, you are using np.radians(10), which a machine float, since that is what NumPy uses. This completely defeats the purpose of computing the final answer to 150 digits, since this input parameter is only accurate to 15. Consider using mpmath.pi/18 instead to get a value that is correct to the number of digits you specified.

user_1618010
Related