Re: Speeding up a numerical integrator
- To: mathgroup at christensen.cybernetics.net
- Subject: [mg1895] Re: Speeding up a numerical integrator
- From: John Tanner <John at janacek.demon.co.uk>
- Date: Sat, 12 Aug 1995 22:52:36 -0400
- Organization: Peace with the World
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.
>
> 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);
>
> .... Note
> that even at 40 the answer is 0.1346, which is still rather far off
> from the correct answer of .125. ......
In fact, the Do function is using the wrong limits and increment. try:
xincr=N[(xmax-xmin)/div]; yincr=N[(ymax-ymin)/div]; zincr=N[(zmax-zmin)/div];
....
{i,xmin+xinxr/2,xmax-xincr/2,xincr},
{j,ymin+yinxr/2,ymax-yincr/2,yincr},
{k,zmin+zinxr/2,zmax-zincr/2,zincr}
using offset increments avoids having to use Simpsons rule etc (orrible).
This now gives answer=0.125 even for div=2 !-))
Have fun,
John.
--
I hate this 'orrible computer : I really ought to sell it
It never does what I want : but only what I tell it. -- (cookie)
John Tanner John at janacek.demon.co.uk 100344.3241 at compuserve.com