Possible bug in numerical integration

[Peter Luschny]

Hi,

I tested on SMC, SageMath 6.10.

Consider the numerical integral:

def T(v):

    def f(t): return (tanh(exp(i*t))/exp(i*t*v)).real()

    c = integral_numerical(f(t), 0, 2*pi)[0]

    return (c*gamma(v+1)/(2*pi)).n()

print [round(T(n)) for n in range(10)]

Sage returned: [0, 1, 0, -1, 0,  8, 0, -136, 0, 3968]

I expected:    [0, 1, 0, -2, 0, 16, 0, -272, 0, 7936]

Looks like there is a factor of 2 missing somewhere.

Checked with Maple:

T := proc(v) local f,c;

f := t -> Re(tanh(exp(I*t))/exp(I*t*v));

c := int(f(t), t=0..2*Pi);

round(evalf(c*GAMMA(v+1)/(2*Pi))) end:

seq(T(n), n=0..9);

Peter

[Peter Luschny]

Further investigation shows that things work with the 

integrand defined as:

def f(t): return exp(x*exp(i*t)-v*i*t)*(exp(exp(i*t))/cosh(exp(i*t))-1)

So the source of the trouble seems to be related to the identity 

   (exp(exp(I*t))/cosh(exp(I*t))-1) = tanh(exp(I*t))

Peter

[Robert Dodier]

On 2016-05-31, Peter Luschny <peter....@gmail.com> wrote:

> def T(v):
> def f(t): return (tanh(exp(i*t))/exp(i*t*v)).real()
> c = integral_numerical(f(t), 0, 2*pi)[0]
> return (c*gamma(v+1)/(2*pi)).n()
>
> print [round(T(n)) for n in range(10)]
>
> Sage returned: [0, 1, 0, -1, 0, 8, 0, -136, 0, 3968]
> I expected: [0, 1, 0, -2, 0, 16, 0, -272, 0, 7936]

With Maxima built from recent source (approximately Maxima 5.38),
I get results that agree with what you expected:

for v:0 thru 9 do
print(float(gamma(v+1)/(2*%pi)*first(quad_qags(f(t),t,0,2*%pi))));
=>
- 4.417437057588218E-17
1.0
- 5.300924469105862E-17
- 1.99999999999998
- 6.361109362927033E-16
15.99999999999981
3.180554681463517E-15
- 271.9999999999945
6.055776113506536E-12
7936.000000000022

My first guess is that the real part was being computed incorrectly and
now it's fixed. For realpart(tanh(exp(%i*t))/exp(%i*t*v)) I get:

(sin(5*t)*sin(2*sin(t)))/(cos(2*sin(t))+cosh(2*cos(t)))
+(cos(5*t)*sinh(2*cos(t)))/(cos(2*sin(t))+cosh(2*cos(t)))

Are you getting something different for that? If so maybe that explains
the difference.

I guess I'm also assuming that Sage punts to Maxima for real() here. But
integral_numerical is probably not calling Maxima, right?

best

Robert Dodier

[Peter Luschny]

> With Maxima built from recent source (approximately Maxima 5.38),
> I get results that agree with what you expected

Thanks Robert, yes, I suspected the numerical integration
prematurely. A plot clearly shows that the culprit are the
real parts of the hyperbolic functions.

plot([tanh(exp(i*t)).real(), (exp(exp(i*t))/cosh(exp(i*t))-1).real()],t,0,2*pi)
The two functions are identical, the plot shows different functions.

> For realpart(tanh(exp(%i*t))/exp(%i*t*v)) I get:
> [...]

> Are you getting something different for that?

I get for (tanh(exp(i*t))/exp(i*t*v)).real()

-e^(imag_part(v)*real_part(t) + imag_part(t)*real_part(v))*
sin(imag_part(t)*imag_part(v) - real_part(t)*real_part(v))*
tan(e^(-imag_part(t))*sin(real_part(t)))/(tan(e^(-imag_part(t))*
sin(real_part(t)))^2*tanh(cos(real_part(t))*e^(-imag_part(t)))^2+1)+
cos(imag_part(t)*imag_part(v) - real_part(t)*real_part(v))*
e^(imag_part(v)*real_part(t) + imag_part(t)*real_part(v))*
tanh(cos(real_part(t))*e^(-imag_part(t)))/(tan(e^(-imag_part(t))*
sin(real_part(t)))^2*tanh(cos(real_part(t))*e^(-imag_part(t)))^2 + 1)

Cheers,
Peter

[kcrisman]

Thanks Robert, yes, I suspected the numerical integration
prematurely. A plot clearly shows that the culprit are the
real parts of the hyperbolic functions.

plot([tanh(exp(i*t)).real(), (exp(exp(i*t))/cosh(exp(i*t))-1).real()],t,0,2*pi)
The two functions are identical, the plot shows different functions.

> For realpart(tanh(exp(%i*t))/exp(%i*t*v)) I get:
> [...]
> Are you getting something different for that?

I get for (tanh(exp(i*t))/exp(i*t*v)).real()

-e^(imag_part(v)*real_part(t) + imag_part(t)*real_part(v))*
sin(imag_part(t)*imag_part(v) - real_part(t)*real_part(v))*
tan(e^(-imag_part(t))*sin(real_part(t)))/(tan(e^(-imag_part(t))*
sin(real_part(t)))^2*tanh(cos(real_part(t))*e^(-imag_part(t)))^2+1)+
cos(imag_part(t)*imag_part(v) - real_part(t)*real_part(v))*
e^(imag_part(v)*real_part(t) + imag_part(t)*real_part(v))*
tanh(cos(real_part(t))*e^(-imag_part(t)))/(tan(e^(-imag_part(t))*
sin(real_part(t)))^2*tanh(cos(real_part(t))*e^(-imag_part(t)))^2 + 1)

This might be Pynac-related.  Numerical integration with that command might be GSL, not Maxima (you may want to try the .nintegrate() method for Maxima). 

[Ralf Stephan]

On Friday, June 3, 2016 at 9:45:02 AM UTC+2, Peter Luschny wrote:

plot([tanh(exp(i*t)).real(), (exp(exp(i*t))/cosh(exp(i*t))-1).real()],t,0,2*pi)
The two functions are identical, the plot shows different functions.

Your Sage is too old, this Pynac bug (existing for years) was fixed

months ago and should be in 7.2.

[Peter Luschny]

> Your Sage is too old, this Pynac bug (existing for years) was fixed
> months ago and should be in 7.2.

Thanks Ralf!

As I said I work on SMC. So I have to wait until William updates the cloud.

Peter

[Ralf Stephan]

In an SMC terminal session:

~$ sage

┌────────────────────────────────────────────────────────────────────┐

│ SageMath Version 6.10, Release Date: 2015-12-18                    │

│ Enhanced for SageMathCloud.                                        │

└────────────────────────────────────────────────────────────────────┘

Really half a year behind, William?

[Dima Pasechnik]

on SMC terminal

$ sage-develop

gives you 7.3.beta2

[Peter Luschny]

> on SMC terminal
> $ sage-develop
>
> gives you 7.3.beta2

Now how can I use this kernel in my yupyter notebook?
Ahh, it's there also! Is this new? Anyway thanks for this hint!

Peter

[Peter Luschny]

This was perhaps not a good idea. It's been long since I got so many ugly error messages.

/projects/sage/sage-dev/src/sage/libs/ecl.pyx in sage.libs.ecl.ecl_safe_eval (/projects/sage/sage-dev/src/build/cythonized/sage/libs/ecl.c:4779)()

    340     if ecl_nvalues > 1:

    341         s = si_coerce_to_base_string(ecl_values(1))

--> 342         raise RuntimeError("ECL says: "+ecl_base_string_pointer_safe(s))

    343     else:

    344         return ecl_values(0)

RuntimeError: ECL says: Detected access to an invalid or protected memory address.

 

[William Stein]

Absolutely. It's far more important that Sage be stable for people
using SMC, especially during the academic semester, than use the
latest version. 6 --> 7 also had a lot of major changes that I wanted
to let shake out, and broke a lot of our integration, and may take a
lot of work to fix.

But, as mentioned elsewhere, we also provide sage-develop = latest
version I can get.

When we switch to Docker containers for projects it will be easier to
support stability and bleeding edge simultaneously, but we're not
there today.

William

--
William (http://wstein.org)