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