On Thu, 2010-09-16 at 15:48 -0700, maldun wrote: > Hi all! > > I'm currently found some time to work again on ticket #8321 again, and > have added some numerical integration tests.
Hi Maldun, I think that improvement of numerical integration would be excellent. I've been somewhat disappointed with the numerical integration in sage in the past. Probably most of the functionality is actually in sage or in components of sage, but it isn't exposed too well. (Or else I have always used the wrong functions, or haven't tried in a while.) I have some more thoughts below. > But I found some major problems, that should be discussed. I said > already to Burcin that I will help out with this, but I think there > are some things that need discussion before I start. > > As far as I have seen sage has the following behavior: calculate the > exact integral if possible, numerical eval this, else call _evalf_. > > But this is not necessarily clever! > An example: > sage: integrate(sin(x^2),x,0,pi).n() > --------------------------------------------------------------------------- > TypeError Traceback (most recent call > last) > ... > /home/maldun/sage/sage-4.5.2/local/lib/python2.6/site-packages/sage/ > rings/complex_number.so in > sage.rings.complex_number.ComplexNumber.__float__ (sage/rings/ > complex_number.c:7200)() > > TypeError: Unable to convert -2.22144146907918 + 2.22144146907918*I to > float; use abs() or real_part() as desired > > this error comes from erf trying to evaluate a complex number (which > by the way should possible because it's an anlytic function...) > But one can also run into numerical troubles! See an example which I > posted some time ago: > http://ask.sagemath.org/question/63/problems-with-symbolic-integration-and-then-numerical-evaluating > Normally numerical integration is done completly numeric! For example > mathematica always distinguish between Integrate, and NIntegrate. > N[Integrate[...]] is not the same as NIntegrate! > In sage it seems there are a few different ways to get a numeric approximation to an integral. There are the top level functions integral and numerical_integral (and their synonyms), and then symbolic expressions have a method nintegrate. So if f is a symbolic expression, 1.) numerical_integral(f, 0, 1) 2.) integral(f, x, 0, 1).n() and 3.) f.nintegrate(x, 0, 1) do 3 different things. Maybe there are also other "top level" ways of definite integration that don't require either an import statement or a direct call to mpmath. I think I would like 1 and 3 to do the exact same thing. I don't know if there is a good reason for them to be different. (I do know that I might often be interested in integrating things that aren't symbolic expressions, though.) I think that the behavior of integral().n() that you list above is reasonable, though. If it fails like in the first example you have above, then that is probably a _different_ bug. And the second example is a type of problem we will probably always have to deal with. (I have a similar example that came up in the "real world" that I might try to dig up because I think it will have the same bad behavior.) A different behavior that I would find reasonable is for integral().n() to _never_ do numeric integration. If integral() can't get an exact answer, then integral().n() could just give an error and ask the user to call numerical_integral() if that's what's wanted. If integral().n() always did numeric integration, then it would of course suffer from the problems below. > Then on the other hand if we would just hand over things to libraries > like pari and mpmath we get no hint if our results are trustworthy. > Consider following simple example of an oscillating integral: > > sage: mpmath.call(mpmath.quad,lambda a: > mpmath.exp(a)*mpmath.sin(1000*a),[0,pi]) > -1.458868938634994559655 > sage: integrate(sin(1000*x)*exp(x),x,0,pi) > -1000/1000001*e^pi + 1000/1000001 > sage: integrate(sin(1000*x)*exp(x),x,0,pi).n() > -0.0221406704921088 > > Do you see the problems?! These are caused by the high oscillation, > but we get no warning. > If you use scipy you would get the following: > > sage: import scipy > sage: import scipy.integrate > sage: scipy.integrate.quad(lambda a: exp(a)*sin(1000*a),0,pi) > Warning: The integral is probably divergent, or slowly convergent. > (-0.1806104043408597, 0.092211010882734368) > > That's how it should be: The user has to be informed if the result > could be wrong, and is then able to choose another method. But in the > current state, bad bad things could happen... > Printed warnings like this can be useful, but they can also be annoying, especially if the numeric integration routine is not called from interactive code. And they might not be necessary in an example like above where the error estimate returned is almost as big as the answer returned. At the very least, if a warning might be printed to standard output I would like a way to turn it off. mpmath's approach might be good, where it seems the error is only returned if requested. Perhaps the default might be to print a warning if the error is not requested and the answer is suspected to be bad, but to not print any warning if the error is going to be returned. > I already posted this issue on the mpmath devel-group. I hope Fredrik > Johansson reads this also. > > Any thoughts or hints about that? > My main other thought is: Thanks for trying to improve numerical integration. Since I've already said so much I may try to look at the code/tickets soon to see if I have an opinion on how things should be done, instead of just opinions on what I would like. > greez, > maldun > -- To post to this group, send an email to sage-devel@googlegroups.com To unsubscribe from this group, send an email to sage-devel+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sage-devel URL: http://www.sagemath.org
