integration not working

[sagetrac-mariah]

sage: g = abs(x^2-1)^(-2/3)
sage: for i in range(3,15):
....:     RR(g.integrate(x, -10^i, 10^i))
....:     
12.0193735393553
12.3400570107857
12.4880961536322
12.5548875258538
12.6007931754768
12.5769884999380
12.6035401600491
12.4014006867096
11.1318879504315
7.82165950876205
6.48361691536824
6.32209292067043
sage: 

Values should increase as g is a positive function!

CC: @kcrisman @burcin

Component: calculus

Reviewer: Karl-Dieter Crisman, Burcin Erocal

Issue created by migration from https://trac.sagemath.org/ticket/10550

[sagetrac-mhampton]

comment:1

Seems like there might be two issues here. One is that the default for integrate_numerical (which is what I think gets called here after symbolic integration fails) only uses a maximum of 87 sample points. This parameter can be set by using integrate_numerical directly but I don't think it can be through integrate.

The second possible issue is that it seems that integrate_numerical tries to create a fast version of the function, but this is not done in an accurate way. I believe this code was written prior to fast_callable, which seems to help:

g = abs(x^2-1)^(-2/3)
gc = fast_callable(g,vars=[x],domain=RDF)
for i in range(3,15):
    print i, RR(g.integrate(x, -10^i, 10^i)), integral_numerical(gc, -10^i, 10^i, max_points=150)

3 12.0193735393553 (12.019608153321734, 9.1957321661164813e-05)
4 12.3400570107857 (12.341098980932561, 9.1150919425152691e-05)
5 12.4880961536322 (12.49028406429874, 6.9889256990153531e-05)
6 12.5548875258538 (12.55959625938481, 7.1073622858611871e-05)
7 12.6007931754768 (12.591761943525762, 0.00010294503436049848)
8 12.5769884999380 (12.60668258943978, 0.00010062716727896715)
9 12.6035401600491 (12.613617891834338, 0.00014081135169109776)
10 12.4014006867096 (12.616777202843952, 9.2533059420607869e-05)
11 11.1318879504315 (12.618201503159618, 0.0025695049437327948)
12 7.82165950876205 (12.618006422456004, 0.00021237998608546343)
13 6.48361691536824 (12.618918276607321, 0.0032425594623566767)
14 6.32209292067043 (12.619569344159935, 0.0097174091356236336)

[sagetrac-mhampton]

comment:2

It might also be worth noting that integrating over the entire real line is special-cased by integrate_numerical, so we get:

RR(g.integrate(x, -Infinity, Infinity))
12.6178509639583

This is a little low, Mathematica gives 12.6195 but its hard to improve the precision. This is an intrinsically hard sort of problem to automate (but we can clearly improve the current behavior).

[williamstein]

comment:3

The reason this is broken is because the implementation of the method _evalf_ in sage/symbolic/integration/integral.py in DefiniteIntegral completely ignores the parent and the precision of a and b. This function has to be rewritten to take into account the precision of the inputs:

    def _evalf_(self, f, x, a, b, parent=None):
        from sage.gsl.integration import numerical_integral
        return numerical_integral(f, a, b)[0]

[williamstein]

comment:4

My personal opinion is that it should work like this:

sage: RR(g.integrate(x, -10^15, 10^15))
error message: you have to specify some params to numerical_integral when constructing the definite integral
sage: RR(g.integrate(x, -10^15, 10^15, max_points=10^5))
correct answer

Or maybe not. I hate this.

[williamstein]

comment:5

Mike Hansen: One option -- just rewrite _evalf_ to call mpmath, which is not perfect, but way better than GSL.

[williamstein]

comment:6

Be sure to use fast_callable:

sage: f(x) = cos(x)^10 - exp(x*sin(cos(x)))
sage: g = fast_callable(f, RealField(200), [x])
sage: timeit('f(10.0)')
^C625 loops, best of 3: 1.25 ms per loop
sage: a = RealField(200)(10)
sage: f(a)
0.17239044387825963285967105367392659309547067375921972300438
sage: g(a)
0.17239044387825963285967105367392659309547067375921972300438
sage: timeit('f(a)')
625 loops, best of 3: 1.38 ms per loop
sage: timeit('g(a)')
625 loops, best of 3: 40.7 µs per loop
sage: 1380/40.7
33.9066339066339

[sagetrac-gagansekhon]

comment:7

some more information about integration in mpmath and other programs

http://groups.google.com/group/mpmath/browse_thread/thread/2fb6d36501c49a73

[sagetrac-mhampton]

comment:8

Your timings with fast_callable seem misleading to me since you are not including the overhead of constructing the fast callable. I got:

sage: f(x) = cos(x)^10 - exp(x*sin(cos(x)))
sage: g = fast_callable(f, RealField(200), [x])
sage: def gmaker(afunc):
    return fast_callable(afunc, RealField(200), [x])
sage: timeit('f(a)')
125 loops, best of 3: 1.46 ms per loop

sage: timeit('g(a)')
625 loops, best of 3: 65.2 µs per loop

sage: timeit('gmaker(f)(a)')
125 loops, best of 3: 1.55 ms per loop

[sagetrac-gbe]

comment:9

Not entirely misleading though, if you call the fast_callable many many times, as would seem to be the case here.

[sagetrac-gbe]

comment:10

Getting this right will be tricky. While simply changing from a blind application of gsl's numerical_integral() to mpmath's quad() may be often be an improvement, in this situation it makes things worse:

sage: f = abs(x^2-1)^(-2/3)
sage: g = fast_callable(f, RealField(100), [x])
sage: import mpmath as mp
sage: mp.mp.prec = 100

sage: for i in [mp.quad(g, [-10^i, 10^1]) for i in range(15)]:
....:     print i
....:     
+inf
8.8396932208278888892057527426
8.3969760821447628003659891605
9.7953936385023669571121865442
9.2339519516756815076612878988
8.2441701346282922281629165563
8.044969834509411578891250552
11.635099287297832304011545535
11.084646393563409097851890886
7.5468622876138167706495612921
7.0790386862087011969090472793
8.4443867135724817987750671703
10.690050592136777734927819129
32.336415464969432639295368434
10.533897172238729354640752465

sage: mp.quad(g, [mp.mpf('-inf'),mp.mpf('+inf')])
mpf('+inf')

Switching from tanh-sinh quadrature to Gauss-Legendre gives

sage: for i in [mp.quadgl(g, [-10^i, 10^1]) for i in range(15)]:
....:     print i
....:     
6.3780650727028379121031555789
7.9564100202151404330674115745
8.3052587750738550159021915092
11.156320092824073496435254839
9.5630214580799381652111853848
2.8857506804320681789970931663
1.7298436143374586397937735747
0.63762636717628812549878116562
0.29045600775700733200460854022
0.13456948161705169202354793393
0.062450129901992433601218328054
0.02898624937362903007167714756
0.013454200394425034174020875196
0.0062448854829133373771884652354
0.002898619019151035832588794256

sage: mp.quadgl(g, [mp.mpf('-inf'),mp.mpf('+inf')])
mpf('11.350773575585495163385287559838')

Hand-tuning with Gauss-Legendre:

sage: mp.quadgl(g, [mp.mpf('-inf'),-1,1,mp.mpf('+inf')])
mpf('12.350413467045610425723223632513')

sage: mp.quadgl(g, [mp.mpf('-inf'),-1,1,mp.mpf('+inf')], maxdegree=12)
mpf('12.592883474959440301710672862944')

And with tanh-sinh:

sage: mp.quad(g, [mp.mpf('-inf'),-1,1,mp.mpf('+inf')])
mpf('+inf')
sage: mp.quad(g, [mp.mpf('-inf'),-1,1,mp.mpf('+inf')], maxdegree=12)
mpf('+inf')

It seems this integral requires at least some hand tuning to be reasonable.

[sagetrac-gbe]

comment:11

There's probably a whole class of functions -- positive functions with a singularities at a finite number of points that go to zero sufficiently fast at plus/minus infinity maybe -- which exhibit this sort of breakage.

[sagetrac-gagansekhon]

comment:12

I tried gp

%gp

for(n=1,10, print(intnum(x=-10^n,10^n,(abs(x^2-1))^(-2/3))))

and got

7.9618235972792581897233298717819451083
6.7165473168753097702622726677655297805
25.355067401739496118829772887176008798
235.26460811139471211432045489286370351
2344.1629174317876121014814723385308094
23437.691659284423942658077351512842284
234375.08896035939931707271806395149682
2343750.0412917410819973945084611231908
23437500.019165928432661951605079616090
234375000.00889603593988653387093897736

[sagetrac-dsm]

comment:13

Just for the record, if given enough precision and the singularity locations, mpmath does okay in tanh-sinh on the infinite range:

sage: import mpmath
sage: dpss = [int(10*1.25**i) for i in [0..15]]
sage: qq = []
sage: for dps in dpss:
....:         mpmath.mp.dps = dps
....:     q = mpmath.quad(lambda x: abs(x^2 - 1)^(-2/3), (-infinity, -1, 1, infinity))
....:     delta = float(abs((q-qq[-1])/qq[-1])) if qq else 0.0
....:     print dps, RealField(250)(q), delta
....:     qq.append(q)
....:  
10 12.619415730000000000000000000000000000000000000000000000000000000000000000 0.000000000000000
12 12.619588240700000000000000000000000000000000000000000000000000000000000000 1.36705461222e-05
15 12.619632936725900000000000000000000000000000000000000000000000000000000000 3.54179740682e-06
19 12.619638678800019440000000000000000000000000000000000000000000000000000000 4.55011184915e-07
24 12.619638942183136436066100000000000000000000000000000000000000000000000000 2.08708920835e-08
30 12.619638947875526360028104307100000000000000000000000000000000000000000000 4.51073913449e-10
38 12.619638947928984872554788554326908474000000000000000000000000000000000000 4.2361364495e-12
47 12.619638947929088225923816071123105700944331096000000000000000000000000000 8.18988320141e-15
59 12.619638947929088350562971317084699893895128992807003256090000000000000000 9.87660231487e-18
74 12.619638947929088350575171596623717955924858573775002931329627429009895074 9.66769302146e-22
93 12.619638947929088350575171711452592280507575074285006026048789429945225951 9.0992202549e-27
116 12.619638947929088350575171711452647219166086056903447254987222698078233522 4.35342554075e-33
145 12.619638947929088350575171711452647219167199848815619049997783218559507282 8.82586195031e-41
181 12.619638947929088350575171711452647219167199848815874865637945653424795957 2.02712328949e-50
227 12.619638947929088350575171711452647219167199848815874865637945883686562894 1.82463038671e-62
284 12.619638947929088350575171711452647219167199848815874865637945883686562894 9.97427124086e-78

and plotting it shows a very plausible logarithmic decrease in the size of the change (can't really call it an error..) You have to give it way more digits of precision than it actually gets right, but it's possible that it converged in the end.

[kcrisman]

comment:15

Burcin, I think that this might be a dupe of #8321, at least for practical purposes. What do you think?

[burcin]

Reviewer: Burcin Erocal

[burcin]

Author: Karl-Dieter Crisman

[burcin]

comment:16

Yep, this is a duplicate. I suggest we close this. I'll write a comment on #8321 so the discussion here is not lost. We were looking for examples on that ticket. It would be great if the examples here could be transferred somehow. Any volunteers?

[jdemeyer]

Changed author from Karl-Dieter Crisman to none

[jdemeyer]

Changed reviewer from Burcin Erocal to Karl-Dieter Crisman, Burcin Erocal