Improvement of numerical integration routines
maldun
I'm currently found some time to work again on ticket #8321 again, and
have added some numerical integration tests.
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!
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...
I already posted this issue on the mpmath devel-group. I hope Fredrik
Johansson reads this also.
Any thoughts or hints about that?
greez,
maldun
Fredrik Johansson
> 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:
It is possible to get an error estimate back from mpmath, as follows:
sage: mpmath.quad(lambda a: mpmath.exp(a)*mpmath.sin(1000*a),
[0,mpmath.pi], error=True)
(mpf('-0.71920642950258207'), mpf('1.0'))
sage: mpmath.quad(lambda a: mpmath.exp(a)*mpmath.sin(1000*a),
[0,mpmath.pi], maxdegree=15, error=True)
(mpf('-0.022140670492108779'), mpf('1.0e-37'))
Currently it doesn't seem to work with mpmath.call (because
mpmath_to_sage doesn't know what to do with a tuple).
The error estimates are also somewhat bogus (in both examples above --
the first value has been capped to 1.0, if I remember correctly what
the code does, and the second is an extrapolate of the quadrature
error alone and clearly doesn't account for the arithmetic error). But
it should be sufficient to correctly detect a problem and signal a
warning in the Sage wrapper code in typical cases.
I'm currently working on rewriting the integration code in mpmath to
handle error estimates more rigorously and flexibly. This will work
much better in a future version of mpmath.
Fredrik
maldun
On 17 Sep., 22:45, Fredrik Johansson <fredrik.johans...@gmail.com>
wrote:
That would be great if this would be fixed!
Don't hesitate to ask me for help. I'm still a beginner, but can
always help out with basic programming
and testing!
maldun
rjf
selected points cannot provide
a rigorous error bound. I don't know what a rigorous "error estimate"
is, unless it is an
estimate of what might be a rigorous bound. And thus not rigorous at
all, since the
estimate might be wrong.
A nasty function: f(x) := if member(x,evaluation_point_set) then 0
else 1.
integrate(f,a,b) should be b-a. The numerical integration program
returns 0.
There are SO many programs to do numerical quadrature, and some of
them are
quite good and free. Some are not free (e.g. NIntegrate in
Mathematica).
The real issue is not writing the program, but figuring out which
program to use and
how to set up the parameters to make it work best. I doubt that
providing
arbitrary-precision floats as a foundation will help solve all the
issues, though
a version of Gaussian quadrature or Clenshaw-Curtis using bigfloats
would
not be hard to do. See http://www.cs.berkeley.edu/~fateman/papers/quad.pdf
for program listings of a few quadrature routines using maxima and
bigfloats.
See what NIntegrate does, and think about
whether you want to do that.
RJF
On Sep 17, 1:45 pm, Fredrik Johansson <fredrik.johans...@gmail.com>
wrote:
maldun
On Sep 19, 5:12 pm, rjf <fate...@gmail.com> wrote:
> Any program that computes a quadrature by evaluating the integral at
> selected points cannot provide
> a rigorous error bound. I don't know what a rigorous "error estimate"
> is, unless it is an
> estimate of what might be a rigorous bound. And thus not rigorous at
> all, since the
> estimate might be wrong.
>
estimation
especially for smooth functions.
But it is the wrong way to say: Hey I can't do it
for every case, so why bother?
The example I gave is smooth simple function but causes already
troubles.
> A nasty function: f(x) := if member(x,evaluation_point_set) then 0
> else 1.
> integrate(f,a,b) should be b-a. The numerical integration program
> returns 0.
that is necessary to allow the user
to give customized integration points (I can show you much less nasty
functions), and for example scipy allows this, and mpmath allows
different integration methods. Of course for special points
error estimation is hard or even impossible, but if the user do it by
him/herself one could assume that he/she knows
what it does.
...
> RJF
>
greez,
maldun
rjf
On Sep 20, 7:14 am, maldun <dom...@gmx.net> wrote:
>
> That's true, but it is important that automated routines do good error
> estimation
> especially for smooth functions.
a program
to test if a function is smooth :)
> But it is the wrong way to say: Hey I can't do it
> for every case, so why bother?
> The example I gave is smooth simple function but causes already
> troubles.
known kind
of problem that can be done -- if you recognize it as oscillatory --
by methods
devised by Filon, or more recently Levin, Iserles, Olver, ... But is
difficult
by regular quadrature programs.
There is a great deal written about integration. Perhaps instead of
offering
to do some basic programming you could offer to read some books or
papers on
numerical integration. Probably the documentation for NIntegrate
would be a start.
maldun
On 21 Sep., 16:23, rjf <fate...@gmail.com> wrote:
> On Sep 20, 7:14 am, maldun <dom...@gmx.net> wrote:
>
>
>
> > That's true, but it is important that automated routines do good error
> > estimation
> > especially for smooth functions.
>
> That is pretty easy if you have a smooth function. So perhaps we need
> a program
> to test if a function is smooth :)
>
> > But it is the wrong way to say: Hey I can't do it
> > for every case, so why bother?
> > The example I gave is smooth simple function but causes already
> > troubles.
>
> The example you gave is what is called "oscillatory" and is a well
> known kind
> of problem that can be done -- if you recognize it as oscillatory --
> by methods
> devised by Filon, or more recently Levin, Iserles, Olver, ... But is
> difficult
> by regular quadrature programs.
>
the comments above from
Fredrik Johansson it seems to me that it was not a such bad test...
Perhaps I should sort out my point before we cause misunderstandings:
It's true that a user familiar with numerics knows about such
behavior. But since
Sage is a program for a wide range of users, I think it's important to
offer user friendly
integration routines also for people who are beginners, or not so
familiar with that topic.
Perhaps I'm wrong, but isn't it a good thing to open sage to a broader
audience?
And of course all packages are included in Sage, so it's no big deal
to use them if I want for example
specific integration points.
And I don't claim here it's some job for overnight. I'm well aware of
this, but I think it should be considered
as long time goal, and I don't see anything wrong about that, and I
believe it deserves some discussion.
And since I offer my own manpower and knowledge to
work on this improvements, I personally think it's nothing wrong about
making suggestions, and point out problems.
But it's quite disturbing for me that something simple like
integrate(sin(x^2),x,0,pi).n() already crashes.
But I think there are already some little improvements, which are not
such big deal, e.g. offer the possibility for the user
to set a specific algorithm in the .n(), or distinguish between
integral.n() and n(integral)
> There is a great deal written about integration. Perhaps instead of
> offering
> to do some basic programming you could offer to read some books or
> papers on
> numerical integration. Probably the documentation for NIntegrate
> would be a start.
not a pro in Python, and will hopefully change that.
I'm just claiming that I can't help out with too advanced programming
skills.
rjf
On Sep 22, 4:12 pm, maldun <dom...@gmx.net> wrote:
>
> Perhaps I should sort out my point before we cause misunderstandings:
> It's true that a user familiar with numerics knows about such
> behavior.
is oscillatory, and it is somewhat of a specialized area to know
about this kind of quadrature.
>But since
> Sage is a program for a wide range of users, I think it's important to
> offer user friendly
> integration routines also for people who are beginners, or not so
> familiar with that topic.
one from a user who is basically ignorant. Giving right answers
to easy questions but wrong answers to difficult questions does
not fit the bill.
> Perhaps I'm wrong, but isn't it a good thing to open sage to a broader
> audience?
> And of course all packages are included in Sage, so it's no big deal
> to use them if I want for example
> specific integration points.
a route for that "broader audience" to use them, or they will not
be used.
If you cannot "integrate" such programs smoothly into a grand
unified program that does numerical integration in Sage, it is
not of much use.
> And I don't claim here it's some job for overnight. I'm well aware of
> this, but I think it should be considered
> as long time goal, and I don't see anything wrong about that, and I
> believe it deserves some discussion.
> And since I offer my own manpower and knowledge to
> work on this improvements, I personally think it's nothing wrong about
> making suggestions, and point out problems.
you are not (I suspect) skilled in numerical integration technology.
So what is wrong with making suggestions and pointing out
problems is that your suggestions may be naive, and the problems
you point out are not actual problems but merely difficulties
YOU have because you are unfamiliar with the standard
solution techniques. So that is what is potentially wrong.
> But it's quite disturbing for me that something simple like
> integrate(sin(x^2),x,0,pi).n() already crashes.
If Maxima is being called, then Sage already has access to the answer.
If Sage screws this up, maybe it has nothing to do with integration.
It has to do with Maxima <--> Sage communication.
> But I think there are already some little improvements, which are not
> such big deal, e.g. offer the possibility for the user
> to set a specific algorithm in the .n(), or distinguish between
> integral.n() and n(integral)
There are substantial numbers of quadrature routines and options
in Maxima, if Sage just allows the user to pass these specs through.
perhaps you should look into what Maxima offers.
>
> > There is a great deal written about integration. Perhaps instead of
> > offering
> > to do some basic programming you could offer to read some books or
> > papers on
> > numerical integration. Probably the documentation for NIntegrate
> > would be a start.
>
> I think to read documentation and papers is a natural thing. But I'm
> not a pro in Python,
and will hopefully change that.
Python
is generally a bad idea.
> I'm just claiming that I can't help out with too advanced programming
> skills.
numerical methods. If you have neither, you should acquire them.
Jonathan Bober
> 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.
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
>
Jonathan Bober
> Any program that computes a quadrature by evaluating the integral at
> selected points cannot provide
> a rigorous error bound. I don't know what a rigorous "error estimate"
> is, unless it is an
> estimate of what might be a rigorous bound. And thus not rigorous at
> all, since the
> estimate might be wrong.
This is, of course, true, provided that no information is known about
the function except for its values on a set of measure 0. However, it is
certainly possible for a numeric integration routine to provide rigorous
error bounds provided that some additional constraints are satisfied.
For example, many functions have bounded variation, or have a bounded
derivative, or are strictly increasing. And I have found situations
where I would like to provide that information to an integration routine
and have it use it in order to provide me information about the error
bound. And it probably isn't unreasonable for an integration routine to
guess at that information, and to provide error bounds that are rigorous
assuming a set of not unreasonable hypotheses about the function it is
integrating.
>
> A nasty function: f(x) := if member(x,evaluation_point_set) then 0
> else 1.
> integrate(f,a,b) should be b-a. The numerical integration program
> returns 0.
>
I am aware that this type of problem can actually occur in real world
examples, but I suspect that in the specific case of numeric integration
it can be somewhat alleviated by using random sampling to choose an
initial evaluation_point_set. Regardless, the existence of unsolvable
problems should not stop us from solving the solvable problems.
Burcin Erocal
Jonathan Bober <jwb...@gmail.com> wrote:
> 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.)
AFAIK, the reason 1 and 3 work differently is mostly historic. Fixing
this could be a part of #7763:
http://trac.sagemath.org/sage_trac/ticket/7763
The behavior of integral().n() has to be different from
numerical_integral() since integral() first tries to integrate
symbolically.
One can use integral(..., hold=True).n() to get a similar effect, but
the options you can pass to the numerical evaluation routines of
symbolic expressions are limited at the moment.
> 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.
I would be OK with this as well.
> If integral().n() always did numeric integration, then it would of
> course suffer from the problems below.
How does doing N[Integrate[ ... ]] in MMA compare to using NIntegrate?
Cheers,
Burcin
Jason Grout
> At the very least, if a warning might be printed to standard
> output I would like a way to turn it off.
If you use the standard Python warning framework, this would be easy.
Jason
rjf
On Sep 22, 11:57 pm, Jonathan Bober <jwbo...@gmail.com> wrote:
....
> guess at that information,
> and to provide error bounds that are rigorous
> assuming
>a set of not unreasonable hypotheses about the function it is
> integrating.
solved.
For example, you could solve World Hunger if you assume
(a) there is enough free food for everyone and (b) there is
enough free transportation to distribute it.
I have no objection to your idea of coming up with an estimate. Just
drop the word "rigorous" .
RJF
rjf
On Sep 23, 5:46 am, Burcin Erocal <bur...@erocal.org> wrote:
> How does doing N[Integrate[ ... ]] in MMA compare to using NIntegrate?
Integrate[f[x],{x,a,b}] exactly by symbolic methods. Then
evaluates the result if possible. NIntegrate does not try to do a
symbolic integral first.
maldun
is to do a design document, or at least start a ticket where we come
to a common agreement how to handle the input because there are
several possiblities to do that. And it is important that we are more
or less happy about the solution.
I'm a little bit short on time in the next weeks, but I hope to get
back into it soon =)
@rjf Of course from a mathematical point of view you are correct, but
from a more pragmatic point of view you have to consider that about
80%-90% of the input satisfy these assumptions. So one can prevent
about 80%-90% of the mistakes. For the others it is important to give
the user a good control about the choice of the method.
It remembers me somehow of the arguments people used against
seatbelts, because there are SOME cases where the seatbelt could kill
you even when it saves your life in 80% of the cases...
>For example, you could solve World Hunger if you assume
>(a) there is enough free food for everyone and (b) there is
>enough free transportation to distribute it.
rjf
On Sep 28, 4:02 am, maldun <dom...@gmx.net> wrote:
> Ok first I would thank everyone for their thoughts. I think the best
> is to do a design document, or at least start a ticket where we come
> to a common agreement how to handle the input because there are
> several possiblities to do that. And it is important that we are more
> or less happy about the solution.
> I'm a little bit short on time in the next weeks, but I hope to get
> back into it soon =)
>
> @rjf Of course from a mathematical point of view you are correct, but
> from a more pragmatic point of view you have to consider that about
> 80%-90% of the input satisfy these assumptions. So one can prevent
> about 80%-90% of the mistakes. For the others it is important to give
> the user a good control about the choice of the method.
program
may be in error by an unknown amount, instead of saying "rigorous
estimate".
Would you be happy if I gave you a cosine routine that gave you a
correct
answer to the last bit except if the input was too close to a rational
multiple of pi, in which
case the error could be anything? Would you argue that a rigorous
estimate
was "correct to the last bit"?
> It remembers me somehow of the arguments people used against
> seatbelts, because there are SOME cases where the seatbelt could kill
> you even when it saves your life in 80% of the cases...
a seatbelt because most trips in a car do not involve a crash.
>
> >For example, you could solve World Hunger if you assume
> >(a) there is enough free food for everyone and (b) there is
> >enough free transportation to distribute it.
>
> Ever watched "We feed the world"? Just curious :)
Dr. David Kirkby
> Ok first I would thank everyone for their thoughts. I think the best
> is to do a design document
I never thought I'd hear anyone used the term "design document" in relation to
Sage!
Seriously, Sage does need to
* Document the requirements for Sage
* Document a specification
* Document how those design requirements will be met
You clearly have hit a problem with numerical integration, but it seems to me
this lack of any design is hurting Sage in many areas. There's no logical path
that gets taken, based on a design.
dave