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)
Change History (19)
comment:1 Changed 3 years ago by dkrenn
- Description modified (diff)
comment:2 Changed 3 years ago by kcrisman
comment:4 Changed 3 years ago by kcrisman
- Description modified (diff)
- Milestone set to sage-4.8
comment:5 follow-up: ↓ 6 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
- Status changed from new to needs_review
comment:8 follow-up: ↓ 9 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: ↓ 10 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!
comment:11 Changed 3 years ago by kcrisman
- 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: ↓ 15 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
A very relevant data point is that
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.
And none of the standard things we do to initialize
change this. And that's not the error message with the documented principal value one that should diverge.
Two more data points of no surprise.
Not sure what is going on here; maybe a missed translation?