Funny behavior with int
289)
My Calculus III students stumbled on this buggy thing while evaluating a line integral to calculate the flux
The curve [X(t),Y(t)] is the right-half of a Lemniscate with polar equation
R^2=cos(2*theta).
The vector field is
F(x,y)=M(x,y)i+N(x,y)j.
They were integrating M*dy-N*dx around the curve.
If we let a=M*dy and aa=expand(M*dy), then they find that Maple's int gives inconsistent results.
As far as I can tell, a and its twin aa are well-behaved over -Pi/4..Pi/4 and equal. Maybe it is a bug in how Maple handles elliptic integrals? Or maybe it is some issue with removable discontinuities?
Code follows: notice how int(a,t=-Pi/4..Pi/4) and
int(aa,t=-Pi/4..Pi/4) are not equal, but ought to be equal.
........
X:=t->cos(t)*sqrt(cos(2*t));
Y:=t->sin(t)*sqrt(cos(2*t));
M:=(x,y)->x^2+4*y;
N:=(x,y)->x+y^2;
a:=M(X(t),Y(t))*diff(Y(t),t);
aa:=expand(a);
int(a,t=-Pi/4..Pi/4);
int(aa,t=-Pi/4..Pi/4);
Int(aa,t=-Pi/4..Pi/4);evalf(%);
plot([a,aa+.1],t=-Pi/4..Pi/4);
- restart:
with(plots):
spacecu
25 min 35 sec ago - to Robert Israel:
yes,
1 hour 59 min ago - to acer:
Yes, it is
2 hours 13 min ago - I agree
2 hours 22 min ago - simpler
2 hours 35 min ago - Not that way
2 hours 53 min ago - Eigenvectors
3 hours 4 min ago - Thanks, acer,
I find another
4 hours 11 min ago - interesting idea
7 hours 14 min ago - interesting
7 hours 17 min ago
- Job offer / Stellenangebot: Technischer Support Maple
- Table from procedure
- Maple Mentor Awards: The Big Kahuna
- The actual architecture
- new package adoption
- Conjectures ripe for Maple experimentation
- overloaded operators
- poll: change in Maple that you want most?
- The future of FourierTrigSeries package
- Maple 12.01 is now Available
Funny behavior with int
I don't know if this is a bug or not. It might have to do with what the function 'expand' is designed to do. I replaced 'expand' with 'normal' and everything turned out as you expected. It could be that your functions evaluate to something with a sqrt(-1) and expand does not handle this condition. See the help pages for 'expand' and 'normal' for a more detailed explanation of this.
Regards,
Georgios Kokovidis
Dräger Medical
bug in int
Let
bbe the first operand ofexpand(a).Maple does not integrate this correctly.
The answer is incorrect, because the integrand is positive everywhere on the open interval. Check:
even more ...
symbolical: [op(aa)]; map( _x -> int(_x,t=-Pi/4..Pi/4), %); `+`(op(%)); evalf(%); [0, 0, 0, 0, 0, 0, 0] 0 numerical: [op(aa)]: map( _x -> evalf(Int(_x,t=-Pi/4..Pi/4)), %); `+`(op(%)); [1.7355011477181, -0.90246059681647, -0.97188064272214, 0.69420045908582, 0., 0., 0.] 0.55536036726531 or Int(a,t=-u..u); value(%); limit(%,u=Pi/4, left); evalf(%); 1/2 2 Pi ------- 16 0.27768018363490But Mathematica gets it correct
Thanks for isolating the problem.
Now I see that even for a simpler situation like
int(cos(t)*sqrt(cos(2*t)),t=-Pi/4..Pi/4)
Maple uses a discontinuous antiderivative, and so gets the definite integral wrong.
Mathematica, on the other hand, returns the continuous antiderivative
(1/4)*arcsin(sqrt(2)*sin(t))*sqrt(2)+(1/2)*sin(t)*sqrt(cos(2*t))
and gets the definite integral correct.
This is kind of depressing.
These issues are now tracked
The integration issues mentioned in this post are now tracked in Maple's internal bug tracking system. Thank you all for identifying the problem and isolating a potential cause for the bug. :)
int
As a long-term solution, perhaps it might be worth considering a full implementation of the Risch algorithm. This came up in another thread. A full Risch implementation would strengthen Maple's integration/DE capabilities; it would also be nice for users to know when an elementary antiderivative did not exist. It would take lots of work, but would be valuable.
Totally useless
For solving problems of this kind, Risch's algorithm is totally useless. In fact, this answer may well be coming out of Risch. The Risch algorithm (in full or in parts) is well-known to not care at all about continuity, and to produce wildly discontinuous answers for clearly continuous inputs. Mathematica is clearly not using (classical) Risch to get this integral correct.
The problem is that the Risch algorithm is purely algebraic, while integration, definite or indefinite, is a problem of analysis. Even though this issue has been identified a long time ago, very few people have made much progress in this area. This is a large gaping hole in the theory, and it does not look like it will be filled any time soon.
My (admittedly radical) opinion: the Fundamental Theorem of Calculus is one of the worst ways of solving problems of definite integration. I used to think FTOC was brilliant as both a theorem and as a computational means; now I have revised that to just 'as a theorem'. This is from years of hard-won experience fixing bugs in definite integration routines. Computationally, computing the limits required to apply FTOC seem to be (in general) a more difficult problem than the original integration problem!
Risch algorithm
Kind thanks for the explanation—I had erroneously understood that things were all solved.
options to FTC?
Jacques,
What other computational approaches are there to definite integration other than the Fundamental Theorem?
I suppose
identify(evalf(Int(f(x),x=a..b));
is one possibility when the result is a number.
The "identify" command is a scary thing, but after this recent experience, I almost feel more comfortable with this than with int(f(x),x=a..b));
You would think that the int routine would do a reality check to see if the result given by FTC via the Risch algorithm makes sense by comparing with evalf(Int...), but I guess it does not.
numerical verification is slow
The problem with numerical verification, done automatically, is that it makes the system much slower.
It's pretty straightforward for the end-user to get a floating-point verification out of Maple, for exact definite integration examples. But if int() always did it, then it'd be measurably slower. Many people might be unhappy about that. Leaving it as a choice for the user, after the computation, is nicer.
The idea is very good for one related purpose, however. It is very good as for assuring that int()'s exact definite integration schemes are behaving and performing well. Produce a scheme for automatically generating "random" exact definite integration test problems, and such a scheme can becomes strong test engine. By "random" I mean selections from a true miscellany.
acer
progress
I wonder, how many of the issues mentioned in Davenport's 2003 Calculemus paper are still relevant today? (It mentions some other methods done by Maple, Alex, such as table-lookup and convolution of MeijerG functions.)
Also, one year after an earlier mapleprimes post, how much of its content is still to the point?
acer
Re: progress
On the theory side, very little progress has been made since that paper [which I quite like, but from the footnotes, you can probably figure out why!]. It is possible that some practical work inside Maplesoft has been done to make things better - I guess we will find out when Maple 12 comes out.
For definite integration, other than table-lookup and methods based on MeijerG functions [both implemented and used heavily], there are also methods based on contour integrals [implemented in Maple but the code is seldom used], integral transforms [not implemented], and differential operators [not implemented].
In the mean time, recheck
While theory and implementation on definite integration progresses, I think that the best practice is not to believe blindly on the results and recheck them in as many ways as possible and practical. Some options are:
1. Use the available symmetries. In this case it should hold
2. Watch the plot of the indefinite integral within the interval of integration for "something strange". In this case the finite discontinuity at the origin.
3. Evaluate the integral with another CAS as different programs/versions seem to use different "recipes". Eg. as pair of additional examples, in Mupad 4:
int(cos(t)*sqrt(cos(2*t)),t);1/4*2^(1/2)*arcsin(2^(1/2)*sin(t)) + 1/2*2^(1/2)*sin(t)*(cos(t)^2 - 1/2)^(1/2)While in Mupad 2:
int(cos(t)*sqrt(cos(2*t)),t);1/2*sin(t)*(1 - 2*sin(t)^2)^(1/2) - 1/4*I*2^(1/2)*ln(I*sin(t)*2^(1/2) + (1 - 2*sin(t)^2)^(1/2))automatic recheck
So if the philosophy is that we should always recheck definite integrals because Maple (like other systems) has a history of incorrectly calculating definite integrals, then it makes sense that Maple should automatically cross-check by comparing with evalf(Int(f(x),x=a..b)) whenever the result is numeric (no left over parameters).
But acer says people would not appreciate this because it would slow down int(f(x),x=a..b). Surely having Maple do the reality check would take less time than having the human remember to follow with evalf(Int(f(x),x=a..b)).
no, please no automatics like that ...
... there are certainly examples, where numerics will fail (large integer exponents or Gamma), where numerics can be terrible expensive ... and IIRC even conversely
check option?
But adding an option to int, at choice of the user, for numerical checking the integral, could perhaps be helpful for a set of integrals at identifying "issues". I wonder whether this conjecture has been explored.
you did have an example some time ago where that would not work
you did have an example some time ago where that would not work, i think it was an integral where liftable singularity (or they have been numerical singular?) have been close to the boundary or similar
edited: ... was it Int(cosh(v)/sqrt(sech(v)-sech(vm)),v=0..vm) ?
Surely you will find examples
But here I am talking about statistics on large number of cases.
If a numerical check were shown to be effective frequently I believe that it could be a useful addition (as everything employed with proper care!)
Unlikely
Such an option, while potentially quite useful [though sometimes misleading as well], has one big problem: it would essentially admit that Maple has bugs. Otherwise, what would be the point of that option?
Of course, if one is serious about applying Maple to problems in engineering, then that option is not a nicety, but manditory. Interesting conundrum, really.
not automatically
It's a good idea to check the final results of any computation, regardless of whether it comes from a CAS or some other program. My emphasis would be on the word "final", there.
A particular definite integral result might well not be the final result. If the software double checks all intermediate results, then it'll be far slower. Why not check only the final result, if it is a compound and more involved computation? If problems are found with the final result, then go ahead and attempt a forced test of intermediate results.
There's another important issue with numerical floating point tests of symbolic results. If you don't know the numerical conditioning of the problem, how are you going to know what fixed precision at which to do the computation? How will you decide, for example, that small imaginary floating-point components are not merely artefacts? Letting the program try to figure this out has at least two major snags. The first is that it may be very difficult to compute the conditioning (and hence also the needed working precision to get a desired accuracy). The second is that if the program is allowed to automatically adjust the precision then it might try to use an absurd number of working digits (arbitrary, since the example is unspecified).
A check of final results does not necessarily have to be a quantitative test, such as a floating-point numeric comparison. It could be a qualitative test, according to some characteristic that is deduced.
All software has bugs. Most all computational software has a history of producing some incorrect results. Even an exact definite integral result, corroborated by a floating-point approximation, should be checked somehow if it's going to be used for some important purpose. Adding an automatic floating-point verificaton attempt, with its resulting performance penalty, should be up to the individual to select. It would be a judgment, weighing the liklihood of error against the performance cost.
acer
example
Digits:=14; # numerical int is ok Int(1/(1+t^2),t=a..b): subs(b=a+1,%): subs(a=5*10^7,%); evalf(%); 50000001 / | 1 | ------ dt | 2 / 1 + t 50000000 -15 0.39999999200000 10 # checking the symbolic does not work int(1/(1+t^2),t=a..b); subs(b=a+1,%); subs(a=5*10^7,%); evalf(%); -arctan(a) + arctan(b) -arctan(a) + arctan(a + 1) -arctan(50000000) + arctan(50000001) 0. # except after increasing digits significantly evalf[40](%%): evalf(%); -15 0.39999999200000 10 # or using an equivalent expression arctan((b-a)/(1+a*b)); subs(b=a+1,%); subs(a=5*10^7,%); evalf(%); b - a arctan(-------) 1 + a b 1 arctan(-------------) 1 + a (a + 1) arctan(1/2500000050000001) -15 0.39999999200000 10 The last comes through the additions theorem i think, that example is taken from: Richard J. Fateman and W. Kahan: Improving Exact Integrals From Symbolic Algebra Systems (2000)Discontinuous Elliptic Substitution
Some of the weirdness seen in int in this thread comes from a subtle issue in the elliptic integrator (and thus it affects a lot of integrals with square-roots of trig functions).
Maple is often pretty good at applying the fundamental theorem of calculus to compute definite integrals even in a case when the indefinite integral it returns is discontinuous.
For the problem above where
f:=2*cos(t)^5*sqrt(2*cos(t)^2-1);a:=-Pi/4; b:=Pi/4;Maple tries a couple different technique. Unfortunately it happens to return the buggy answer:which as I mentioned, comes from an issue with trig substitution in the elliptic integrator:
`int/definite/elliptic`(f,t,a,b,{}); # 0 `int/definite/ftoc`(f,t,a,b,{}); evalf(%); # 1/2 # 25 2 Pi # ---------- # 64 # # 1.735501148As you can see, there is an internal integrator "ftoc" that gets the right answer. However, "ftoc" it is a lot slower than "elliptic" (1/10s vs. 2s on my machine) and that probably gives you a hint as to why "elliptic" is the first algorithm tried here.
Simple question? (maybe)
My, apologies if this seem like a simple question. Why does Maple, or any other CAS, return a discontinuous anti-derivative in situations where a continuous anti-derivative is guaranteed to exist by FTC? My problem is that I can think of some reason why this might be done but I am only guessing since I do not know for certain. (I won't bother listing any since they may be wrong)
Thanks,
Thomas
Constructive algorithms
FTC is a non-constructive theorem - it only tells you that there exists a continuous anti-derivative, not how to find it. Current algorithms for finding an anti-derivative, the famous Risch algorithm and its many extensions, is well-known to ignore issues of continuity altogether, since it is a purely algebraic algorithm. No one has yet designed an algorithm for integration in finite terms that respects analysis, except for some very small families of expressions.
Thanks
Jacques, thanks for the reply. I had not really paid close attention to these issues prior to today. I mostly use Maple to do textbook type problems since I teach HS. Now I know I can't just compute an anti-derivative and take it for granted that I can apply FTC. Strangely, despite the fact that this issue has been discussed before it didn't really sink in until I read this thread in conjunction with jpmay's recent blog A Perilous Tale of Definite Integration.
Simple answer?
"Guaranteed to exist" is not the same thing as "easy to calculate". The continuous antiderivative that the FTC guarantees is the definite integral from some fixed origin to x. But that's not useful if you want a symbolic antiderivative (especially if you'll want to use that symbolic antiderivative to calculate the definite integral).
For example, consider f(x) = max(0, g(x)) where g has a known, continuous antiderivative G. So a (discontinuous) antiderivative of f is piecewise(g(x) > 0, G(x), 0). But to construct a continuous antiderivative of f you'd need to know all the points where g changes sign.
Thanks
Robert, thanks for the reply. I knew it was only an existence theorem, but I was unsure how Maple or other CAS appoach the problem. For example, I was wondering if they first attempted to find a continuous antiderivative. I can see now that that way of thinking is a gross oversimplification of the problem in general.
Fixed in Maple 12
For what it is worth, the issue reported in Alex's original post was fixed in Maple 12:
> int(a,t=-Pi/4..Pi/4); 1/2 2 Pi ------- 8 > int(aa,t=-Pi/4..Pi/4); 1/2 2 Pi ------- 8> evalf(Int(a,t=-Pi/4..Pi/4)); identify(%); 0.5553603673 1/2 2 Pi ------- 8 > evalf(Int(aa,t=-Pi/4..Pi/4)); 0.5553603673John
value
It's worth a lot to know that bugs are tracked and corrected. Thanks.