Is there anyway to reduce the time for the numerical integration of this formula?

Question

I'm trying to do the following series of computations.

data = Import["http://www.vibrationdata.com/elcentro_EW.dat"];
Acc = Interpolation[data];
c = 1500; rho = 1; H = 80;
p[y_, t_] := (4*rho*c)/π*
   Sum[(-1)^(n - 1)/(2 n - 1.) Cos[(2 n - 1.) π/(2 H) y]*
     NIntegrate[
      BesselJ[0, (2 n - 1) π/(2 H)*c*(t - τ)] * Acc[τ], 
       {τ, 0, t}], {n, 1, 20}];

(* the following takes a long time and I would like to make it fast *)
(* Table[{p[10, t]}, {t, 0, 8, .1}]; *)

The list of data is shown below. The data can be found here:http://www.vibrationdata.com/elcentro_EW.dat

And the constants values are: c=1500, rho=1, H=80.

Is there anyway to accelerate this computation? I want to evaluate p at y=10.

Thank you.

Answer 1

To be really fast, you bypass NIntegrate. But you can only do that if you know ahead of time what sort of sampling you need, since you lose two great advantages of NIntegrate, adaptive sampling and error estimation. In this case, the integral consists cubic polynomials times a Bessel function over many, small intervals of the same width. (The cubic polynomials come from the interpolating function Acc). Well, this is not a hopeless situation for Gauss-Legendre integration. The oscillations of the Bessel function in a given interval increase linearly with n in the OP's formula. One expects to have to add a couple of Gauss points for each oscillation. Next one needs a starting number of points for n = 1. I used 5 at first but then I started testing (not shown*) the accuracy and found that 3 is sufficient for 7-8 digits precision.

Clear[p, gaussrule, grdata];
mem : grdata[n_] := mem = NIntegrate`GaussRuleData[n, MachinePrecision];
gaussrule[f_, gp_, intervals_?(MatrixQ[#, NumericQ] &)] := 
  Module[{abscissae, weights},
   {abscissae, weights} = Most@grdata[gp];
   f@Map[
      #[[1]] + abscissae (#[[2]] - #[[1]]) &,
      intervals
      ].weights.Subtract[intervals[[All, 2]], intervals[[All, 1]]]
   ];

c = 1500; rho = 1; H = 80;
p[y_, t_] := 
  With[{intervals = Partition[Append[TakeWhile[data[[All, 1]], # < t &], N@t], 2, 1]},
   (4*rho*c)/π*
    Sum[
     (-1)^(n - 1)/(2 n - 1) Cos[(2 n - 1) π/(2 H) y]*
      gaussrule[
       Function[τ, BesselJ[0, (2 n - 1) π/(2 H)*c*(t - τ)] * Acc[τ]],
       Ceiling[3 + n/2],  (* Gauss points heuristic *)
       intervals],
     {n, 1, 20}]
   ];

I said "really fast" at the beginning, but that's only in comparison with the OP's code:

Table[{p[10, t]}, {t, 0, 8, .1}]; // AbsoluteTiming
(* {29.541, Null}  *)

For comparison, for p[10, 8], this method takes 0.67 s., the OP's method 31.4 s., and Method -> "InterpolationPointsSubdivision" 16.8 s.

*Note: Testing can be done by observing the convergence of the computed values as the number of points is increased, for instance. I did not do a thorough testing. The accuracy is inferred from the theory of convergence of Gauss-Legendre integration and a few spot checks. The 7-8 digits was determined by a couple of high-precision tests using the same approach as gaussrule but with greater number of points and 32-digit precision numbers. These checks took less than ten seconds and were still much faster than NIntegrate.

Another note: BesselJ unpacks packed arrays. Consequently, I didn't try to keep things packed: you often lose some time when they get packed, because they are then unpacked.