MathGroup Archive 1995

[Date Index] [Thread Index] [Author Index]

Search the Archive

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


  • Prev by Date: Re: Functional programming puzzle
  • Next by Date: Re: Mathematica for Linux is out
  • Previous by thread: MS-DOS Version
  • Next by thread: Speeding up a numerical integrator