9
$\begingroup$

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]

Mathematica graphics

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}}
$\endgroup$

2 Answers 2

10
$\begingroup$

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]

plot

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.

$\endgroup$
1
  • $\begingroup$ Thanks a lot. I had tried something similar but my s values were off. Looking at your much better approach now, I see why. $\endgroup$
    – smayhem
    Commented Sep 27, 2014 at 18:55
7
$\begingroup$

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] &

Mathematica graphics

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.