# Efficient way to obtain values of a function defined by an Integral

Consider the following equation:

$$S(q)=\frac{(4 \pi \rho ) \int r (h(r)-1) \sin (q r) \, dr}{q}$$

I want to numerically obtain values for $S(q)$ given that I have data points representing $h(r)$ over the limits of interest of the integral above (in this case $[0.05, 11.4]$ ). Here is what I have so far and it works but seems slow. I first created an interpolation of $h(r)$

h = Interpolation[hrdata]; (* see below for hrdata *)


Then used the following function (with $\rho = 0.83$ ):

s[q_?NumericQ] := 4 Pi 0.83 NIntegrate[r Sin[q r] (h[r] - 1), {r, 0.05, 11.4}] / q


This Table creates the desired output:

ans = Table[{i, s[i]}, {i, 0.05, 11.4, 0.1}];


This is the plot:

ListLinePlot[ans, PlotRange -> All]


So my question is, is there a better more efficient way to do this even if NIntegrate doesn't have to be used?

Here is the data used:

hrdata = {{0.0478682458, 0.}, {0.143604737, 0.}, {0.239341229, 0.}, {0.335077721, 0.},
{0.430814212, 0.}, {0.526550704, 0.}, {0.622287196, 0.}, {0.718023687, 0.},
{0.813760179, 0.}, {0.909496671, 0.}, {1.00523316, 0.}, {1.10096965, 0.},
{1.19670615, 0.}, {1.29244264, 0.}, {1.38817913, 0.}, {1.48391562, 0.},
{1.57965211, 0.}, {1.6753886, 0.}, {1.7711251, 0.}, {1.86686159, 0.},
{1.96259808, 0.}, {2.05833457, 0.}, {2.15407106, 0.}, {2.24980755, 0.},
{2.34554405, 0.}, {2.44128054, 0.}, {2.53701703, 0.}, {2.63275352, 0.},
{2.72849001, 0.}, {2.8242265, 0.}, {2.919963, 0.}, {3.01569949, 0.0000676622354},
{3.11143598, 0.00451295794}, {3.20717247, 0.075857987}, {3.30290896, 0.425311056},
{3.39864545, 1.18839185}, {3.49438195, 2.10924971}, {3.59011844, 2.73332417},
{3.68585493,  2.94676433},  {3.78159142, 2.82338311}, {3.87732791, 2.52318174},
{3.9730644, 2.19565845}, {4.0688009,  1.86881686}, {4.16453739, 1.58405181},
{4.26027388, 1.35851141}, {4.35601037, 1.16444392}, {4.45174686, 1.00457511},
{4.54748335, 0.905797344}, {4.64321985, 0.815426641}, {4.73895634, 0.739436997},
{4.83469283, 0.686695997}, {4.93042932, 0.655832514}, {5.02616581, 0.624168951},
{5.1219023, 0.598239362}, {5.2176388, 0.590977257}, {5.31337529, 0.591387319},
{5.40911178,  0.602061955}, {5.50484827, 0.616496255}, {5.60058476, 0.634917018},
{5.69632125, 0.671428594}, {5.79205775, 0.703768473}, {5.88779424, 0.757543323},
{5.98353073, 0.815278632}, {6.07926722, 0.871993974}, {6.17500371, 0.92784548},
{6.2707402, 0.984303451}, {6.3664767, 1.04019019}, {6.46221319, 1.10732974},
{6.55794968, 1.14081452}, {6.65368617, 1.18276983}, {6.74942266, 1.21864694},
{6.84515915, 1.2463683}, {6.94089565, 1.26732968}, {7.03663214, 1.2714219},
{7.13236863, 1.27410128}, {7.22810512, 1.25860771}, {7.32384161, 1.23118455},
{7.4195781, 1.19937153}, {7.5153146, 1.15440898}, {7.61105109, 1.11170797},
{7.70678758,  1.06629363}, {7.80252407, 1.0167857}, {7.89826056, 0.978978348},
{7.99399705, 0.94260256}, {8.08973355,  0.916606926}, {8.18547004, 0.881495983},
{8.28120653, 0.864364397}, {8.37694302, 0.849789475}, {8.47267951, 0.842009732},
{8.568416, 0.835537063}, {8.6641525, 0.837210014}, {8.75988899, 0.84566796},
{8.85562548,  0.859726517}, {8.95136197, 0.876656885}, {9.04709846, 0.897456184},
{9.14283495, 0.917005224}, {9.23857145, 0.942830462}, {9.33430794, 0.971618284},
{9.43004443, 0.98840516}, {9.52578092, 1.02075017}, {9.62151741, 1.03977905},
{9.7172539, 1.06007292}, {9.8129904, 1.07522679}, {9.90872689, 1.08518106},
{10.0044634, 1.09122059}, {10.1001999, 1.09955184}, {10.1959364, 1.10425954},
{10.2916729, 1.0968197}, {10.3874093, 1.08912877}, {10.4831458, 1.08586879},
{10.5788823, 1.07263336}, {10.6746188, 1.05734122}, {10.7703553, 1.04470208},
{10.8660918, 1.02973102}, {10.9618283, 1.01332368}, {11.0575648, 1.00350336},
{11.1533013,  0.987646017}, {11.2490378, 0.97643616}, {11.3447743,
0.964238364}, {11.4405108, 0.953274271}}


The interpolation step seems to be unnecessary because the integral into which it enters can be equally well approximated as a Riemann sum. So to get really fast results you could do the following:

{r, h} = Transpose[hrdata];

d = Differences[r];

Clear[s];
s[q_] := (4 Pi 0.83 )/q Total[d Rest[r Sin[q r] (h - 1)]]

ans2 = Table[{i, s[i]}, {i, 0.05, 11.4, 0.1}];

plot2 = ListLinePlot[ans2, PlotRange -> All]


Here I defined r and h as single lists, and took d to contain the differences in the r values so I can approximate the integral as a sum of rectangles with base d. The entries in d happen to be practically identical because the data are apparently on a uniform grid, so this step could be simplified further. But with these lists, the original formula for s[q] can now be written down without any explicit looping constructs because all the required operations automatically thread over lists. To deal with the end points, I simply threw out the first data point because it's zero anyway.

To plot the result I collect the points ans2 and plot them as in the question. The plot is pretty much the same as with interpolation.

• Thanks a lot. I had tried something similar but my s values were off. Looking at your much better approach now, I see why. Sep 27, 2014 at 18:55

Here is the same approach as Jens, but using Association and assuming uniform grid data:

st[hr_?MatrixQ, rho_, min_, max_, step_] :=
Module[{h = <|Rule @@@ hr|>, s, freq, dr = hr[[2, 1]] - hr[[1, 1]]},
s = Function[q, 4 Pi rho /q Tr[dr # Sin[q #] (h[#] - 1) & /@ Keys[h]]];
freq = Range[min, max, step];
Transpose[{freq, s /@ freq}]
]


Now:

st[hrdata, 0.83, 0.05, 11.4, 0.05] // ListLinePlot[#, PlotRange -> All] &