# Ways to quicken numeric integral over interpolating spline?

I have an integral of the form $$\int_1^{15000} f_1(k) f_2(k) dk,$$ where $f_1(k)$ is an interpolating function and $f_2(k)$ is some multiplicative function of k. I'm working with the following

   data = Table[{i, AsRunDec[asMz /. NumDef, Mz /. NumDef, i, 2]}, {i,1, 15000, 14999/100000}];

f1 = Interpolation[data];


Because of the large integration domain, a naive use of NIntegrate seems to give rise to a long run time.

   NIntegrate[f1[k]*f2[k],{k,1,15000}];


I think the longevity is due to mathematica looping over the table of values in data over and over again as it performs the integration. I just wondered what amendments I can make to quicken this evaluation? General methods involving the $f_i$'s and what I've defined above are ok or, if preferred I can post my whole code.

I should say RunDec is a specialised package. Basically for each value of $i$ it returns some value.

Thanks!

Edit: Here is the code I'm working with. Basically I'd like to amend it to make the NIntegrate over the interpolating function (here called cubic1) quicker. In the notation above, $f_2(k)$ depends on some unknown parameters NN,a,b which are determined through a minimisation but probably those details are irrelevant below - I just include for completeness.

   Needs["ErrorBarPlots"]

MasData5 = {{{44.8, 47.5}, ErrorBar[4.0]}, {{54.8, 50.1},
ErrorBar[4.2]}, {{64.8, 61.7}, ErrorBar[5.1]}, {{74.8, 64.8},
ErrorBar[5.5]}, {{84.9, 75}, ErrorBar[6.2]}, {{94.9, 81.2},
ErrorBar[6.7]}, {{104.9, 85.3}, ErrorBar[7.1]}, {{119.5, 94.5},
ErrorBar[7.5]}, {{144.1, 101.5}, ErrorBar[8.3]}, {{144.9, 101.9},
ErrorBar[10.9]}, {{162.5, 117.8},
ErrorBar[12.8]}, {{177.3, 130.2},
ErrorBar[13.4]}, {{194.8, 147.7},
ErrorBar[17.1]}, {{219.6, 137.4},
ErrorBar[20.1]}, {{244.8, 176.6},
ErrorBar[20.3]}, {{267.2, 178.7},
ErrorBar[21.1]}, {{292.3, 200.4}, ErrorBar[29.1]}, {{60, 55.8},
ErrorBar[4.838]}, {{80, 66.6}, ErrorBar[7.280]}, {{100, 73.4},
ErrorBar[6.426]}, {{120, 86.7}, ErrorBar[7.245]}, {{140, 104},
ErrorBar[12.083]}, {{160, 110}, ErrorBar[16.279]}, {{42.5, 43.8},
ErrorBar[3.482]}, {{55, 57.2}, ErrorBar[3.980]}, {{65, 62.5},
ErrorBar[4.614]}, {{75, 68.9}, ErrorBar[5.197]}, {{85, 72.1},
ErrorBar[5.523]}, {{100, 81.9}, ErrorBar[5.368]}, {{117.5, 95.7},
ErrorBar[6.277]}, {{132.5, 103.9}, ErrorBar[6.912]}, {{155, 115},
ErrorBar[7.920]}, {{185, 129.1}, ErrorBar[9.192]}, {{215, 141.7},
ErrorBar[10.666]}, {{245, 140.3}, ErrorBar[14.526]}, {{275, 189},
ErrorBar[24.274]}, {{49, 39.2}, ErrorBar[10]}, {{86, 75.7},
ErrorBar[14.414]}, {{167, 118}, ErrorBar[22.828]}, {{43.2, 50.7},
ErrorBar[1.5]}, {{50, 59.5}, ErrorBar[1.4]}, {{57.3, 61.8},
ErrorBar[1.9]}, {{65.3, 67.6}, ErrorBar[1.7]}, {{73.9, 72.4},
ErrorBar[1.9]}, {{83.2, 79.9}, ErrorBar[2.3]}, {{93.3, 84.4},
ErrorBar[2.1]}, {{104.3, 86.7}, ErrorBar[2.7]}, {{47.9, 55.4},
ErrorBar[2.1]}, {{68.4, 66.4}, ErrorBar[2.9]}};
(*h1 2006 Q^2=0 data,zeus 2002,zeus 2004 and h1 2013 data for Q^2=0*)

gamma = 5.55*^-6;
MJpsi = 3.1;
alphaem = 1/137;
lambda = 0.09;
Ca = 3; (*list of constants *)

xg = NN*((4*qbar)/((4*qbar - MJpsi^2) + w^2))^(-a)*(qbar)^b*Exp[Sqrt[
16*Ca/9*Log[(4*qbar - MJpsi^2 + w^2)/(4*qbar)]*
Log[Log[qbar/lambda]/Log[qo/lambda]]]];

data = Table[{i, AsRunDec[asMz /. NumDef, Mz /. NumDef, i, 2]}, {i, 1,
15000, 14999/100000}];

cubic1 = Interpolation[data];

F5[w_?NumericQ, a_?NumericQ, b_?NumericQ, NN_?NumericQ] :=Module[{qbarr = 2.4025, qoo = 2, qbar, qo}, (3.89379*^5*1/(4.9 + 4*0.06*Log[w/90])/16*4*Pi^3*MJpsi^3*
gamma/12/
alphaem*(NIntegrate[
cubic1[k]/qbar/(qbar + k)*
D[(2^(2*(D[Log[xg],
Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]]
+ 3))/Sqrt[Pi]*
Gamma[D[Log[xg],
Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]]
+ 5/2]/
Gamma[D[Log[xg],
Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]]
+ 4])*
NN*(((4*qbar))/((4*qbar - MJpsi^2) + w^2))^(-a)*(k)^b*
Exp[Sqrt[
16*Ca/9*Log[(4*qbar - MJpsi^2 + w^2)/(4*qbar)]*
Log[Log[k/lambda]/Log[qo/lambda]]]], k] /.
qbar -> qbarr /. qo -> qoo, {k,
qoo, (w^2 - MJpsi^2)/4}] +
0.118/qbar/qo*Log[(qbar + qo)/qbar]*
NN*((4*qbar)/((4*qbar - MJpsi^2) + w^2))^(-a)*(qo)^
b*(2^(2*(D[Log[xg],
Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]]
+ 3))/Sqrt[Pi]*
Gamma[D[Log[xg],
Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]]
+ 5/2]/
Gamma[D[Log[xg],
Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]]
+ 4]))^2) /. qbar -> qbarr /. qo -> qoo]

chisq5[a_?NumericQ, b_?NumericQ, NN_?NumericQ] := Sum[((MasData5[[k, 1, 2]] - F5[MasData5[[k, 1, 1]], a, b, NN])/
MasData5[[k, 2, 1]])^2, {k, 1, Length[MasData5]}];

NMinimize[chisq5[a, b, NN], {a, b, NN}];

• You are getting some good answers but it would be best to supply the full details because as indicated (see Michael E2 answer) the answers are specific to the problem and don't apply in general. – Jack LaVigne Nov 25 '17 at 1:17
• @Jack LaVigne Sure, I just uploaded my whole code I'm working with in the OP. I took the lower boundary to be 2 in this case but it doesn't really matter for my purpose, 1 is also ok. Thanks :) – CAF Nov 25 '17 at 9:24
• I copied the code starting with Needs["ErrorBarPlots"] and it doesn't run. To make your questions easy for members to answer, a good technique for quality control is to close Mathematica. Then start Mathematia and and in the new notebook copy and run your code copies from the question you have written. If it fails, find the mistake or omission and fix it. At a minimum the definitions for AsRunDec and NumDef appear to be missing. – Jack LaVigne Nov 25 '17 at 15:41
• @Jack LaVigne ah yes, sorry RunDec is a specialised Mathematica package and a RunDec.m file needs to be downloaded from the net. Then NumDef will also be understood. – CAF Nov 25 '17 at 16:23
• @JackLaVigne It is RunDec and here: ttp.kit.edu/Progdata/ttp00/ttp00-05 Thanks! – CAF Nov 25 '17 at 19:17

My experience teaches me that most problems tend to be specific, not general. Here's a way that works for a similar looking problem:

f1 = Interpolation@Table[{x, Sin[10 Sin[10 x]]}, {x, 0., 1., 2.^-21}];

NIntegrate[f1[k]*Exp[k], {k, 0, 1}] // AbsoluteTiming

(*  {2.77656, 0.0636131}  *)


100 times faster by hand:

npts = 11;
err = 1;
int = 1;
PrintTemporary@Dynamic@{npts, err, int};
(While[
Abs[err/int] > 10^-10,
npts = 2 npts + 1;
{int, err} = Block[{abscissae, weights, errweights},
{abscissae, weights, errweights} =
NIntegrateGaussKronrodRuleData[npts, MachinePrecision];
{weights, errweights}.(f1[k]*Exp[k] /. k -> abscissae)
]
];
int) // AbsoluteTiming

(*  {0.023723, 0.0636131}  *)


It works because the singularities are negligible and the oscillation moderate.

• I think I am able to follow your hand procedure. I think it would be helpful if you could write a short paragraph explaining the strategy of your code. In particular, I am unfamiliar with NIntegrateGaussKronrodRuleData although I can guess what it doesn. When I use 100,000 points rather than 2,000,000 I do not experience an increase in speed. – Jack LaVigne Nov 25 '17 at 0:20
• @JackLaVigne I'm not sure it's worth working on until the OP responds...Your example is not the OPs nor is mine. See this tutorial and this one for the use of the rule data. – Michael E2 Nov 25 '17 at 1:12
• Thank you for the references. Very helpful! – Jack LaVigne Nov 25 '17 at 1:18

I like turning these types of integrals into an ODE. Using @MichaelE2's example:

NDSolveValue[
{int'[k] == f1[k] Exp[k], int[0]==0},
int[1],
{k, 0, 1}
] //AbsoluteTiming


{0.139189, 0.0636131}

You can build a simple Riemann sum, if there are a lot of data points and the curve is smooth enough. This is also very fast.

tab = Table[{x, Sin[10 Sin[10 x]]}, {x, 0., 1., 2.^-21}];

f1 = Interpolation@tab;


Using FullForm[f1], you see, that you can extract the original data points from the InterpolatingFunction. (If you dont know the original data)

xval = First@f1[[3]]; yval = f1[[4, 3]];


Calculate half of the difference of two neighboring x-values.

rot = (Drop[RotateLeft[xval], -1] - Rest@RotateRight[xval])/2;


Calculate the multiplicative function

expxval = Exp[xval];


Build the Riemann sum

Total[{Total[rot*Drop[yval, -1]*Drop[expxval, -1]],
Total[rot*Rest[yval]*Rest[expxval]]}] // AbsoluteTiming

(*     {0.0780848, 0.0636131}    *)
`

This works for equal and for varying x-distances.