Re: Speeding up a numerical integrator
- To: mathgroup at christensen.cybernetics.net
- Subject: [mg1805] Re: Speeding up a numerical integrator
- From: rknapp (Robert Knapp)
- Date: Thu, 17 Aug 1995 00:01:36 -0400
- Organization: The Santa Fe Institute
In article <DD6Gp4.402 at wri.com> Scott.A.Hill at williams.edu (Lancelot) writes:
>
> Hello. I've found it necessary to write a simple numerical
> integrator in Mathematica, as NIntegrate is not cooperating with me
> for some reason. ...
Perhaps you would do well to contact technical support at WRI for
assistance in using NIntegrate for your particular problem. There are
cases NIntegrate cannot handle, but there are techniques which can
help in many other cases.
> ...Anyway, here is the "Integrator"
>
> Integrator::usage="Integrator[f,{x,xmin,xmax},{y,ymin,ymax},{z,zmin,zmax},div]
> does a numerical integration over the variables x, y, and z, to the
> limits specified, with div divisions per dimension.";
>
> Integrator[f_,{x_,xmin_,xmax_},{y_,ymin_,ymax_},{z_,zmin_,zmax_},div_]:=
> (
> answer=0.;
> Do[(
> part=(f/.{x->i,y->j,z->k})(xmax-xmin)(ymax-ymin)(zmax-zmin)/div^3;
> answer=N[answer+If[NumberQ[part],part,0]]
> ),
> {i,xmin,xmax,N[(xmax-xmin)/div]},
> {j,ymin,ymax,N[(ymax-ymin)/div]},
> {k,zmin,zmax,N[(zmax-zmin)/div]}
> ];
> answer);
>
> I have the "If" statement there in case the function blows up
> to infinity at some point. I figured it was faster than using a Limit
> command on everything.
>
This slows the convergence of your method even below first order near
a singularity. If you are interested in accurate results and you
expect to have singularities, this is very unlikely to speed things
up.
> However, the thing runs rather slowly. For the command
> Integrator[x*y*z,{x,0,1},{y,0,1},{z,0,1},20] it took approximately 28
> seconds. Replacing the 20 by 40 upped the time to 205 seconds. Note
> that even at 40 the answer is 0.1346, which is still rather far off
> from the correct answer of .125. At this rate the thing will take
> 4000 seconds = over an hour to evaluate
> Integrator[x*y*z,{x,0,1},{y,0,1},{z,0,1},100], and my function is much
> more complex than x*y*z.
In a numerical program for a given problem, there are typically at
least two ways of speeding it up. First, to improve the algorithm,
and second to improve the efficiency of the programming. In the
function Integrator above, improvements of both types can easily be
made. I will focus on problems with the algorithm because this is
where the worse problems lie.
Even for a smooth function, the convergence to a solution is
asymptotically only first order in the reciprocal of "div". i.e., the
error goes roughly as c/div. Based on your result, the error is
roughly .384/div. Thus, based on your timings, to get, say 6 digits
of accuracy in the solution (the default for NIntegrate), it would
take div = .384/10^(-6) = 384000. Since the time it will take to do
the problem scales as the cube of div, to get this accuracy it would
take about
In[20]:=
(205.*(384000/40)^3)/3600/24/365
Out[20]=
6
5.75123 10
YEARS to compute on your computer!!!!!
Near a singularity, your method has even slower convergence, so
accurate results would take even more generations.
Reimann sums are a fine theoretical tool, but relatively poor for
numerically approximating integrals. Unless your function has no
smoothness whatsoever, there are better methods. One set of these is
contained in NIntegrate. If it is indeed the case that your problem
is too difficult for NIntegrate to handle (with some help from WRI
technical support), then you may have to use a different method, or do
some analysis to figure out the contibutions to your integral near
singularities and use NIntegrate elsewhere..
Rob Knapp
WRI