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

asked 2019-07-24 14:53:06 +0100

99 ●7 ●10 ●16

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)
   1654         if operator == add_vararg:
   1655             operator = _operator.add

/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)
   1654         if operator == add_vararg:
   1655             operator = _operator.add

/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 flag offensive close merge delete

add a comment

2 Answers

Sort by » oldest newest most voted

answered 2019-07-24 23:49:20 +0100

nbruin
4153 ●3 ●45 ●91

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.

edit flag offensive delete link

Comments

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

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

add a comment

answered 2019-07-24 15:49:25 +0100

Emmanuel Charpentier
7774 ●6 ●50 ●153

updated 2019-07-24 16:26:21 +0100

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)

edit flag offensive delete link

Comments

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.

Fusion809 ( 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).

Emmanuel Charpentier ( 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.

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

add a comment

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

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

2 Answers

Comments

Comments

Your Answer

Question Tools

Stats

Related questions

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

2 Answers

Comments

Comments

Your Answer

Question Tools

Stats

Related questions

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