MathGroup Archive 2002

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

Search the Archive

Re: A faster alternative to ListIntegrate?

  • To: mathgroup at smc.vnet.net
  • Subject: [mg36497] Re: A faster alternative to ListIntegrate?
  • From: Selwyn Hollis <slhollis at earthlink.net>
  • Date: Mon, 9 Sep 2002 00:29:43 -0400 (EDT)
  • References: <3D7B95CC.2000504@earthlink.net>
  • Sender: owner-wri-mathgroup at wolfram.com

For what it's worth, here's a version of Simpson's rule that works with 
both equally spaced and unequally spaced points.


simp = Compile[{x1,y1,x2,y2,x3,y3},
With[{h13=x1-x3, h12=x1-x2, h23=x2-x3},
(-y1*h23*(2x1-3x2+x3) - y2*h13^2 + y3*h12*(x1-3x2+2x3))*h13/(6h12*h23)]];
    
ListSimpson[vals_?VectorQ, dx_Real]:= With[{n=Length[vals]},
If[OddQ[n], (dx/3)*(First[vals] + Last[vals] + 
4*Plus@@vals[[Range[2,n-1,2]]] + 2*Plus@@vals[[Range[3,n-2,2]]]), $Failed]];

ListSimpson[datapts_?MatrixQ] := If[OddQ[Length[datapts]],
    Plus@@(simp@@Flatten[#]&/@ Partition[datapts,3,2]), $Failed]
   

(* equally spaced points*)

data = Table[100*Sin[x], {x, 0, 100, 0.001}];
ListSimpson[data, .001] // Timing

    {0.21 Second, 13.7681}

(* randomly spaced points *)

Table[
{data = Sort[Table[With[{x = 
Random[Real,{0,100}]},{x,100*Sin[x]}],{100001}]];
ListSimpson[data] // Timing,
With[{xvals = Transpose[data][[1]]},
gaps = Drop[RotateLeft[xvals] - xvals, -1];
{Min[gaps], Max[gaps]}]},
{4}]

  {{{1.62 Second, 13.8348}, {1.32223*10^-8, 0.0114454}},
   {{1.59 Second, 13.8261}, {2.42163*10^-8, 0.0127803}},
   {{1.62 Second, 13.8455}, {1.08923*10^-9, 0.0119059}},
   {{1.54 Second, 13.7987}, {9.66814*10^-9, 0.0106988}}}


---
Selwyn Hollis


> Mathew,
>
> Some possibilities
>
>      <<NumericalMath`ListIntegrate`
>
>     ListIntegrate[data]//Timing
>
>         {6.59 Second,13.7681}
>
> The following is suggested in the Help Browser entry for the package
>     Integrate[
>         Interpolation[data, InterpolationOrder\[Rule]1][x],
>     {x,0,100}]//Timing
>
>         {4.56 Second,13.768}
>
> Trapezium rule with equal steps:
>
>     #[[1]]+#[[-1]]+ 2 Tr[Take[#,{2,-2}]]&[data[[All,2]]] 0.01/2//Timing
>
>         {0.22 Second,13.768}
>
> Trapezium rule with possibly unequal steps
>
>     (Drop[#1,1] - Drop[#1,-1]).(Drop[#2,-1] + Drop[#2,1])&[
>       data[[All,1]], data[[All,2]]]/2//Timing
>
> {0.83 Second,13.768}
>
> --
> Allan
>
> ---------------------
> Allan Hayes
> Mathematica Training and Consulting
> Leicester UK
> www.haystack.demon.co.uk
> hay at haystack.demon.co.uk
> Voice: +44 (0)116 271 4198
> Fax: +44 (0)870 164 0565
>
>
> "Matthew Rosen" <mrosen at cfa.harvard.edu> wrote in message
> news:ahr122$l2v$1 at smc.vnet.net...
> > Hi Everyone;
> >  I've tracked down the slow operation of my Mathematica simulation 
> code to
> > lie in the ListIntegrate command:
> >
> > G[n_] := ListIntegrate[xsec Phi[n], 1]
> >
> > where both xsec and Phi[n] are 400 values long.
> >
> > Is there a way to speed up ListIntegrate via Compile or a similar 
> technique?
> >
> > Thanks in advance and best regards,
> >
> > Matt
> > ---
> > Matthew Rosen
> > Harvard-Smithsonian Center for Astrophysics
> > Mail Stop 59
> > 60 Garden Street
> > Cambridge, MA 02138
> >
> > e: mrosen at cfa.harvard.edu
> > o: (617) 496-7614
> >
>
>




  • Prev by Date: RE: RE: RE: Generating Two Unit Orthogonal Vectors
  • Next by Date: Re: Re: Generating Two Unit Orthogonal Vectors (correction)
  • Previous by thread: Getting hierarchies of partitions (continue)
  • Next by thread: Plotting with logarithmic scale?