Opened 3 years ago

Closed 3 years ago

#11987 closed defect (fixed)

integrate returns divergent, maxima.integrate the correct result

Reported by: dkrenn Owned by: burcin
Priority: major Milestone: sage-4.8
Component: calculus Keywords: symbolic integration, maxima, integrate, divergent
Cc: nbruin, mjo Merged in: sage-4.8.alpha4
Authors: Nils Bruin, Karl-Dieter Crisman Reviewers: Karl-Dieter Crisman, Nils Bruin
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Description (last modified by kcrisman)

The following integration fails:

sage: integrate(sin(x)*sin(x/3)/x^2, x, 0, oo)       
ERROR: An unexpected error occurred while tokenizing input
The following traceback may be corrupted or invalid
The error message is: ('EOF in multi-line statement', (554, 0))

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
...
ValueError: Integral is divergent.

But the integral is not divergent. It seems that integrate uses Maxima, but Maxima can compute the correct result:

sage: maxima.integrate(sin(x)*sin(x/3)/x^2, x, 0, oo)
%pi/6

Apply trac_11987-principal-value.patch and trac_11987-doctest.patch.

Attachments (2)

trac_11987-principal-value.patch (1.2 KB) - added by nbruin 3 years ago.
trac_11987-doctest.patch (908 bytes) - added by kcrisman 3 years ago.
Doctests, apply on top of other patch

Download all attachments as: .zip

Change History (19)

comment:1 Changed 3 years ago by dkrenn

  • Description modified (diff)

comment:2 Changed 3 years ago by kcrisman

A very relevant data point is that

sage: maxima_calculus.integrate(sin(x)*sin(x/3)/x^2, x, 0, oo)
    624                 self._name = parent._create(value, name=name)
    625             except (TypeError, KeyboardInterrupt, RuntimeError, ValueError), x:
--> 626                 raise TypeError, x
    627 
    628     def _latex_(self):

TypeError: ECL says: Divergent Integral

And that is the "copy" of Maxima we use for integration and other calculus stuff, and that's why the error is returned.

So the problem must be in one of the "extra" things we use for calculus.


Here's vanilla Sage Maxima.

(%i1) integrate(sin(x)*sin(x/3)/x^2, x, 0,inf);
                                      %pi
(%o1)                                 ---
                                       6

And none of the standard things we do to initialize

init_code = ['display2d : false', 'domain : complex', 'keepfloat : true',
            'load(to_poly_solver)', 'load(simplify_sum)']

change this. And that's not the error message with the documented principal value one that should diverge.

(%i12) integrate(1/x^3,x,-1,3);
Principal Value
(%o12) 4/9

Two more data points of no surprise.

sage: maxima_eval(([max_integrate],[sr_to_max(SR(a)) for a in [sin(x)*sin(x/3)/x^2, x, 0,oo]]))
---------------------------------------------------------------------------
RuntimeError: ECL says: Divergent Integral
sage: ([max_integrate],[sr_to_max(SR(a)) for a in [sin(x)*sin(x/3)/x^2, x, 0,oo]])
([<ECL: $INTEGRATE>], [<ECL: ((MTIMES) ((MEXPT) $X -2) ((%SIN) ((MTIMES) $X ((RAT) 1 3))) ((%SIN) $X))>, <ECL: $X>, <ECL: 0>, <ECL: $INF>])

Not sure what is going on here; maybe a missed translation?

comment:3 Changed 3 years ago by jdemeyer

  • Milestone sage-4.7.3 deleted

Milestone sage-4.7.3 deleted

comment:4 Changed 3 years ago by kcrisman

  • Description modified (diff)
  • Milestone set to sage-4.8

comment:5 follow-up: Changed 3 years ago by kcrisman

  • Cc nbruin added

Nils, do you think that line 88 in maxima_lib.py could conceivably be responsible for these?

ecl_eval('(defun principal nil (error "Divergent Integral"))')

I don't really understand what's going on with that line. Otherwise I don't see how we could be getting something different, though.

comment:6 in reply to: ↑ 5 Changed 3 years ago by nbruin

Replying to kcrisman:

Nils, do you think that line 88 in maxima_lib.py could conceivably be responsible for these?

ecl_eval('(defun principal nil (error "Divergent Integral"))')

You can replace it by this to get closer to the behaviour "throw an error if Maxima wants to print "Principal Value" and do not alter the behaviour otherwise:

ecl_eval('(defun principal nil (cond ($noprincipal (diverg)) ((not pcprntd) (merror "Divergent Integral"))))')

It fixes the example in this report.

However, when you go and take a look in defint.lisp there are some very dodgy (setq pcprintd t) statements in there. Essentially what they are doing there is suppressing the printing of the "Principal Value" message. I hope those pieces of code have established that taking a principal value integral is indeed undeniably the thing to do ...

I think it's defensible to simply patch maxima_lib with this line. We're depending on maxima for integration anyway, so if their handling leads to errors, we can just point "upstream".

For reference defint.lisp, line 1096.

Changed 3 years ago by nbruin

comment:7 Changed 3 years ago by nbruin

  • Authors set to Nils Bruin
  • Status changed from new to needs_review

comment:8 follow-up: Changed 3 years ago by kcrisman

  • Description modified (diff)
  • Reviewers set to Karl-Dieter Crisman
  • Status changed from needs_review to needs_work
  • Work issues set to doctests

Well, needs work of course, as we want a doctest that this works.

I think you are right that this is ok to do. At least I was right that this is what was happening. So what happens, are we replacing their handling of principal with ours in certain cases? I wouldn't mind a quick tutorial in what your function (defun defines one, right?) actually does :)


Various uninformed remarks...

So is this even really a principal value situation? I mean, the 1/x^2 component and then taking each pos. or neg. piece can make the integral from 0 to oo be like an alternating series, in effect... For what it's worth, W|A doesn't even mention PV, just returns pi/6.

There is some interesting numerical instability in integration when you go really far out, not surprising I guess with that crazily oscillating of an integrand.

sage: for end in [2,4,6,8,10,12,14,16]:
    numerical_integral(sin(x)*sin(x/3)/x^2,0,10^end,max_points=10000)
....:     
(0.52351589061571202, 9.9920072216264089e-15)
(0.5235987756295446, 8.2436845585377863e-07)
(0.52359878933232351, 9.998230889899859e-07)
(0.52359885121913463, 1.259313851474051e-06)
(0.52359877780154374, 1.2600489386016923e-06)
(1.2189490243409529e-08, 2.420112684966894e-08)
(7.0214579963551362e-11, 1.4528575097897933e-10)
(1.1777735784870202e-12, 2.382542114112869e-12)

comment:9 in reply to: ↑ 8 ; follow-up: Changed 3 years ago by nbruin

Replying to kcrisman:

So what happens, are we replacing their handling of principal with ours in certain cases?

See the link to defint.lisp above. You can see the original definition of principal there:

1096 (defun principal nil

definition of a function with an empty argument list

1097   (cond

this is the main conditional in lisp. It's like a switch statement

($noprincipal (diverg))

If the variable $noprincipal (the $ is a name-mangling convention in Maxima that this is the maxima variable noprincipal) has a non-nil value, then call diverg and return its value (it actually doesn't return but throws an exception or raises an error)

1098         ((not pcprntd)

If pcprntd has a nil (=false) value:

1099          (format t "Principal Value~%")

Print a message

1100          (setq pcprntd t))))

and set the variable to t to ensure that the message only gets printed once per integral evaluation.

I wouldn't mind a quick tutorial in what your function (defun defines one, right?) actually does :)

All that the new version does is replacing the format print statement by an error-raising statement. We don't set pcprntd because that would not have any effect anyway.

The nasty bits happen in line 1144 and line 1157 where pcprntd gets set to t (note the dynamical scoping here! This is different from what would happen in python. Outside the let statement, pcprntd reverts to its old value again).

For some reason, Maxima decides to evaluate the integral in this report via a principal value integral somewhere, but apparently the writer of the code is confident that this is the right thing to do and suppresses the printing of the Principal Value warning.

With the patch, we accept that if Maxima wasn't going to print the message, things should be safe and no error should be raised.

So is this even really a principal value situation? I mean, the 1/x^2 component and then taking each pos. or neg. piece can make the integral from 0 to oo be like an alternating series, in effect... For what it's worth, W|A doesn't even mention PV, just returns pi/6.

The fact that the patch makes a difference shows that maxima *does* do something that triggers calling principal and that only happens in line 625 and line 1168. The latter is the code path we end up in. Maxima doesn't print the message for the reason explained above.

Go ahead and add any appropriate doctests. Then you can co-author the patch! [and thanks for shaming me into repairing this :-) I hope we don't have other monkey-patch surprises]

comment:10 in reply to: ↑ 9 Changed 3 years ago by kcrisman

Thanks for your explanations, very helpful!

For some reason, Maxima decides to evaluate the integral in this report via a principal value integral somewhere, but apparently the writer of the code is confident that this is the right thing to do and suppresses the printing of the Principal Value warning.

With the patch, we accept that if Maxima wasn't going to print the message, things should be safe and no error should be raised.

That makes sense, actually. The PV value is the value if the integral is defined, as I understand it, so they are using trickery :) Since we aren't going to write our own integration code anyway, we should trust Maxima and report upstream if it does have an error (which this one doesn't).

Go ahead and add any appropriate doctests. Then you can co-author the patch! [and thanks for shaming me into repairing this :-) I hope we don't have other monkey-patch surprises]

Good, I'll do this as I have opportunity over the next couple days. Thanks for figuring this out!

Changed 3 years ago by kcrisman

Doctests, apply on top of other patch

comment:11 Changed 3 years ago by kcrisman

  • Authors changed from Nils Bruin to Nils Bruin, Karl-Dieter Crisman
  • Description modified (diff)
  • Status changed from needs_work to needs_review

Ok, done. Needs review. Passes tests on the symbolic, functions, and calculus directories for me.


Apply trac_11987-principal-value.patch and trac_11987-doctest.patch.

comment:12 Changed 3 years ago by kcrisman

Ping. Should be very easy to finish review, just need the doctest patch to be reviewed.

comment:13 follow-up: Changed 3 years ago by nbruin

  • Status changed from needs_review to positive_review

I think KDC thinks the patch itself is OK and I think his doctest is good and informative, so I think we can give this a positive review between the two of us. Anybody in favour of more red tape should feel free to revert this, though.

comment:14 Changed 3 years ago by kcrisman

  • Reviewers changed from Karl-Dieter Crisman to Karl-Dieter Crisman, Nils Bruin

comment:15 in reply to: ↑ 13 Changed 3 years ago by mjo

  • Cc mjo added

Replying to nbruin:

I think KDC thinks the patch itself is OK and I think his doctest is good and informative, so I think we can give this a positive review between the two of us. Anybody in favour of more red tape should feel free to revert this, though.

I'm not yet qualified to understand the original patch, but I too have reviewed the doctest. I confirmed the result in Mathematica, and ran the long test suite after applying both patches.

comment:16 Changed 3 years ago by jdemeyer

  • Work issues doctests deleted

comment:17 Changed 3 years ago by jdemeyer

  • Merged in set to sage-4.8.alpha4
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.