Complex numerical integration
Oscar Gerardo Lazo Arjona
I've been trying to solve this integral:
sage: integrate(sqrt(sec(x)-1),x,pi/2,pi)
integrate(sqrt(sec(x) - 1), x, 1/2*pi, pi)
Since sage failed to produce a symbolic result, I tried with
numerical_integral
sage: numerical_integral(sqrt(sec(x)-1),pi/2,pi)
(nan, nan)
But that failed also... So I tried with mathematica:
sage: mathematica_console()
Mathematica 7.0 for Linux x86 (32-bit)
Copyright 1988-2008 Wolfram Research, Inc.
In[1]:= Integrate[Sqrt[Sec[x]-1],{x,Pi/2,Pi}]
Out[1]= I Pi
So mathematica says it's a pure imaginary number, which I know to be true.
I know that symbolic integration is quite complicated to develop, but
numerical integration should be fairly easy to extend to complex
numbers. Am I missing something? Why does numerical_integral return NaN?
Perhaps i'ts not because of the complex result, in which case how could
I solve the integral without recurring to mathematica?
thanks!
Oscar.
Johan Grönqvist
> I've been trying to solve this integral:
>
> (nan, nan)
>
> But that failed also... So I tried with mathematica:
>
> numbers. Am I missing something?
A workaround seems to be to integrate the real and imaginary parts
separately:
sage: numerical_integral(real(sqrt(sec(x)-1)),pi/2, pi)
(1.9175999157365625e-16, 5.0010185963949996e-17)
sage: numerical_integral(imag(sqrt(sec(x)-1)),pi/2, pi)
(3.1415926269162875, 2.3498999460392589e-06)
... giving the same result as Mathematica (although numerically).
It seems to me that the integration is performed by gsl, and the gsl
integration routine seems only to take real-valued functions as
arguments, but I have not looked at this previously, so I may be missing
something as well.
Regards
Johan
mhampton
sage: from mpmath import *
sage: mp.dps = 15; mp.pretty = True
sage: f = lambda x: sqrt(sec(x)-1)
sage: quad(f, [pi/2, pi])
(1.43051518370573e-8 + 3.14159264808335j)
To get back to a sage type you could do:
sage: ans=quad(f, [pi/2, pi])
sage: CC([ans.real, ans.imag])
-Marshall
On Oct 13, 11:01 pm, Oscar Gerardo Lazo Arjona
maldun
other thread are are several other problems with numerical integration
which occour.
The thing is, that sage holds a lot of tool to perform numerical
integration to get a good answer (pari, mpmath, scipy etc.) but
sometimes it would be nice to get something more user friendly.
I'm planning to do some work on this, at least improving more control
of the choice of numerical integration methods for the user.
But it won't happening in the next time because I have a lot of work
to do, and I have to finish some other things, like a new version of
orthogonal polynomials before.
Another goal would be doing something similar to Mathematica's
NIntegrate, but his would take some years and should only considered
as long time goal for sage.
-maldun
On Oct 14, 6:01 am, Oscar Gerardo Lazo Arjona
Oscar Lazo
On Oct 14, 4:54 am, Johan Grönqvist <johan.gronqv...@gmail.com> wrote:
> A workaround seems to be to integrate the real and imaginary parts
> separately:
>
> sage: numerical_integral(real(sqrt(sec(x)-1)),pi/2, pi)
> (1.9175999157365625e-16, 5.0010185963949996e-17)
> sage: numerical_integral(imag(sqrt(sec(x)-1)),pi/2, pi)
> (3.1415926269162875, 2.3498999460392589e-06)
the likes of:
sage: numerical_integral(sqrt(sec(x)-1),pi/2, pi,complex=True)
so that it would do the following internally:
rea=numerical_integral(rea(sqrt(sec(x)-1)),pi/2, pi)
ima=numerical_integral(imag(sqrt(sec(x)-1)),pi/2, pi)
return ( rea[0]+I*ima[0], rea[1], ima[1])
what do you think? I could do this myself very easily if we agree that
this is a good idea.
Oscar
Robert Bradshaw
+1, I've had to do that manually myself before.
On Oct 14, 2010 2:40 PM, "Oscar Lazo" <estadist...@gmail.com> wrote:
On Oct 14, 4:54 am, Johan Grönqvist <johan.gronqv...@gmail.com> wrote:
> A workaround seems to be ...
I think it would be enough to add an option to numerical_integral on
the likes of:
sage: numerical_integral(sqrt(sec(x)-1),pi/2, pi,complex=True)
so that it would do the following internally:
rea=numerical_integral(rea(sqrt(sec(x)-1)),pi/2, pi)
ima=numerical_integral(imag(sqrt(sec(x)-1)),pi/2, pi)
return ( rea[0]+I*ima[0], rea[1], ima[1])
what do you think? I could do this myself very easily if we agree that
this is a good idea.
Oscar
--
To post to this group, send an email to sage-...@googlegroups.com
To unsubscribe from this gr...
David Kirkby
IMHO, it would be sensible to make that the default method, which is
what happens in Mathematica's NIntegrate command. You don't need to
specify Nintegrate to return both the real and imaginary parts - it
does that automatically. In this case, you can see it gets a small
real part too, which you know is just rounding errors.
In[1]:= NIntegrate[Sqrt[Sec[x]-1],{x,Pi/2,Pi}]
Out[1]= 3.45315 10^-10 + 3.14159 I
Dave
David Kirkby
> IMHO, it would be sensible to make that the default method, which is
> what happens in Mathematica's NIntegrate command. You don't need to
> specify Nintegrate to return both the real and imaginary parts - it
> does that automatically. In this case, you can see it gets a small
> real part too, which you know is just rounding errors.
>
> In[1]:= NIntegrate[Sqrt[Sec[x]-1],{x,Pi/2,Pi}]
>
> Out[1]= 3.45315 10^-10 + 3.14159 I
>
> Dave
>
On second thoughts, this is not "rounding errors" but likely to be a
combination of both rounding errors and errors introduced as a result
of a fact one breaks the integral into a finite number of sections.
It's probably the latter effect which dominates the actual rounding
errors.
Dave
rjf
No, it doesn't make sense because, in general, separating the real
and the imaginary parts, symbolically
rjf
if what you are doing is separating the real and imaginary parts
symbolically,
then it will not always work at the boundaries of numerical
evaluation.
Consider that you may think that expression xyz is non-negative and
thus
sqrt(xyz) is real. But numerical evaluation makes xyz very very
slightly
negative.
If numerical_integrate(), or rea() does not evaluate its arguments in
the normal fashion,
then that might change the situation.
RJF