# Complex numerical integration: how can it be performed in SageMath?

z=1; numerical_integral(exp(i*(x*z+(x^3)/3)), 0, 1, max_points=100) returns the error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-3-a00a879c8f90> in <module>()
----> 1 z=Integer(1); numerical_integral(exp(i*(x*z+(x**Integer(3))/Integer(3))), Integer(0), Integer(1), max_points=Integer(100))

/usr/lib/python2.7/site-packages/sage/calculus/integration.pyx in sage.calculus.integration.numerical_integral (build/cythonized/sage/calculus/integration.c:2943)()
279                 to_sub = dict(zip(vars[1:], params))
280                 func = func.subs(to_sub)
--> 281             func = func._fast_float_(str(vars[0]))
282         except (AttributeError):
283             pass

/usr/lib/python2.7/site-packages/sage/symbolic/expression.pyx in sage.symbolic.expression.Expression._fast_float_ (build/cythonized/sage/symbolic/expression.cpp:62032)()
11928         """
11929         from sage.symbolic.expression_conversions import fast_float
> 11930         return fast_float(self, *vars)
11931
11932     def _fast_callable_(self, etb):

/usr/lib/python2.7/site-packages/sage/symbolic/expression_conversions.pyc in fast_float(ex, *vars)
1704         1.4142135623730951
1705     """
-> 1706     return FastFloatConverter(ex, *vars)()
1707
1708 #################

/usr/lib/python2.7/site-packages/sage/symbolic/expression_conversions.pyc in __call__(self, ex)
224             return self.tuple(ex)
225         else:
--> 226             return self.composition(ex, operator)
227
228     def get_fake_div(self, ex):

/usr/lib/python2.7/site-packages/sage/symbolic/expression_conversions.pyc in composition(self, ex, operator)
1679         """
1680         f = operator
-> 1681         g = [self(_) for _ in ex.operands()]
1682         try:
1683             return f(*g)

/usr/lib/python2.7/site-packages/sage/symbolic/expression_conversions.pyc in __call__(self, ex)
216                 div = self.get_fake_div(ex)
217                 return self.arithmetic(div, div.operator())
--> 218             return self.arithmetic(ex, operator)
219         elif operator in relation_operators:
220             return self.relation(ex, operator)

/usr/lib/python2.7/site-packages/sage/symbolic/expression_conversions.pyc in arithmetic(self, ex, operator)
1651             from sage.functions.all import sqrt
1652             return sqrt(self(operands[0]))
-> 1653         fops = map(self, operands)

/usr/lib/python2.7/site-packages/sage/symbolic/expression_conversions.pyc in __call__(self, ex)
215             if getattr(self, 'use_fake_div', False) and (operator is _operator.mul or operator is mul_vararg):
216                 div = self.get_fake_div(ex)
--> 217                 return self.arithmetic(div, div.operator())
218             return self.arithmetic(ex, operator)
219         elif operator in relation_operators:

/usr/lib/python2.7/site-packages/sage/symbolic/expression_conversions.pyc in arithmetic(self, ex, operator)
1651             from sage.functions.all import sqrt
1652             return sqrt(self(operands[0]))
-> 1653         fops = map(self, operands)

/usr/lib/python2.7/site-packages/sage/symbolic/expression_conversions.pyc in __call__(self, ex)
206         except TypeError as err:
207             if 'self must be a numeric expression' not in err.args:
--> 208                 raise err
209
210         operator = ex.operator()

TypeError: unable to coerce to a real number


Likewise:

z=1; integral(exp(i*(x*z+(x^3/3))), x, 0, +Infinity)


returns:

integrate(e^(1/3*I*x^3 + I*x), x, 0, +Infinity)


. If I add a .n() after it, to numerically evaluate it, I get much the same error as I showed earlier. So, how can I do numerical integration of complex integrands in SageMath? I know that simple complex integrands (e.g. exp(i*x)) can be analytically integrated by SageMath, but numerical integration does not seem to be an option.

edit retag close merge delete

Sort by ยป oldest newest most voted

I think you can use something along the lines of

CC=ComplexBallField(100)
f=fast_callable(e^(I*(x+x^3/3)),vars=[x],domain=CC)
CC.integral(lambda x,_: f(x)  ,CC(0),CC(1))


which gets you a result along the lines of:

[0.7778906934510079480835322919 +/- 3.85e-29] + [0.5105422937539419671472249270 +/- 1.82e-29]*I


Exciting part: these error bounds are meant to be certified! No other computer algebra package will give you certified numerical integrals. This is using ARB, see http://fredrikj.net/blog/2017/11/new-...

Caveat: this is rather new code, so the interface hasn't been tested very much yet. In particular, while the code above should propagate intervals properly through f, it could easily be the case that some funny intermediate coercion step loses a ball radius somewhere. The integrator does depend on correct propagation of ball radii.

more

Hey, thanks for the link on ARB. I've started reading about it, and I'm totally geeking out!

( 2020-07-21 21:19:11 +0100 )edit

you may try:

fr(x)=(e^(i*x*z+x^3/3)).subs(z==1).real_part()
fi(x)=(e^(i*x*z+x^3/3)).subs(z==1).imag_part()
ir=numerical_integral(fr(x),0,1)
ii=numerical_integral(fi(x),0,1)
Res=tuple([ir[0]+I*ii[0], ir[1]+I*ii[1]])


Would Your Grace accept Res=(0.904174612747833 + 0.5253056923680366*I, 1.0038354733541787e-14 + 5.832064746336536e-15*I) as an approximate answer ?

*Later : Wups, wrong function (forgot a pair of parentheses).

fr(x)=(e^(i*(x*z+x^3/3))).subs(z==1).real_part()
fi(x)=(e^(i*(x*z+x^3/3))).subs(z==1).imag_part()
ir=numerical_integral(fr(x),0,1)
ii=numerical_integral(fi(x),0,1)
Res=tuple([ir[0]+I*ii[0], ir[1]+I*ii[1]])
Res
(0.777890693451008 + 0.5105422937539419*I, 8.636321585109385e-15 + 5.668158095705665e-15*I)

more

Well, it's a roundabout way of solving it, but sure it's satisfactory. Although, if anyone has a more direct way of numerical integration with complex numbers I'll be happy to accept such an answer.

( 2019-07-24 15:52:17 +0100 )edit

As far as I can tell, there's no way to do that in Sage symbolically. Some parts of Sage may be able to do this numerically, but I'm currently not aware of them.

Interestingly, Mathematica cannot do that symbolycally, but can do it numerically (via NIntegrate).

( 2019-07-24 16:16:26 +0100 )edit

This specific integral isn't symbolically integratable at all; I originally meant to write the bounds as -infinity to infinity, in that case the solution is the Airy $2\pi Ai(z)$ function. Complex integrals that are analytically integratable (e.g. $\int_{0}^{2\pi} e^{ix} dx$) can be analytically (or symbolically) solved in SageMath, but not numerically solved.

( 2019-07-24 16:21:35 +0100 )edit