6
$\begingroup$

I have a simple integration which, when using an interpolation function, is taking too long to calculate:

c = 2.99792*10^5;
A = 3.87624*10^-14;
FreeElectronFractionData = {{3000, 1.0829044}, {2984.9246, 1.0828562}, {2969.8493, 1.0828473}, {2954.7739, 1.0828366}, {2939.6985, 1.0828238}, {2924.6231, 1.0828083}, {2909.5478, 1.0827898}, {2894.4724, 1.0827674}, 
    {2879.397, 1.0827404}, {2864.3217, 1.0827077}, {2849.2463, 1.0826683}, {2834.1709, 1.0826207}, {2819.0955, 1.0825632}, {2804.0202, 1.0824939}, {2788.9448, 1.0824106}, {2773.8694, 1.0823111}, {2758.7941, 1.0821927}, 
    {2743.7187, 1.0820531}, {2728.6433, 1.08189}, {2713.5679, 1.0817016}, {2698.4926, 1.0814865}, {2683.4172, 1.0812441}, {2668.3418, 1.0809745}, {2653.2664, 1.0806783}, {2638.1911, 1.0803569}, {2623.1157, 1.0800119}, 
    {2608.0403, 1.0796454}, {2592.965, 1.0792594}, {2577.8896, 1.0788561}, {2562.8142, 1.0784377}, {2547.7388, 1.0780061}, {2532.6635, 1.0775631}, {2517.5881, 1.0771101}, {2502.5127, 1.0766486}, {2487.4374, 1.0761797}, 
    {2472.362, 1.0757042}, {2457.2866, 1.0752228}, {2442.2112, 1.074736}, {2427.1359, 1.0742441}, {2412.0605, 1.0737472}, {2396.9851, 1.0732455}, {2381.9098, 1.0727387}, {2366.8344, 1.0722267}, {2351.759, 1.0717092}, 
    {2336.6836, 1.0711857}, {2321.6083, 1.070656}, {2306.5329, 1.0701194}, {2291.4575, 1.0695754}, {2276.3822, 1.0690234}, {2261.3068, 1.0684627}, {2246.2314, 1.0678928}, {2231.156, 1.0673129}, {2216.0807, 1.0667222}, 
    {2201.0053, 1.06612}, {2185.9299, 1.0655055}, {2170.8545, 1.064878}, {2155.7792, 1.0642365}, {2140.7038, 1.0635802}, {2125.6284, 1.0629083}, {2110.5531, 1.0622197}, {2095.4777, 1.0615136}, {2080.4023, 1.060789}, 
    {2065.3269, 1.0600449}, {2050.2516, 1.0592803}, {2035.1762, 1.0584941}, {2020.1008, 1.0576853}, {2005.0255, 1.0568526}, {1989.9501, 1.0559951}, {1974.8747, 1.0551114}, {1959.7993, 1.0542004}, {1944.724, 1.0532609}, 
    {1929.6486, 1.0522915}, {1914.5732, 1.051291}, {1899.4979, 1.0502581}, {1884.4225, 1.0491914}, {1869.3471, 1.0480894}, {1854.2717, 1.0469509}, {1839.1964, 1.0457743}, {1824.121, 1.044558}, {1809.0456, 1.0433006}, 
    {1793.9703, 1.0420003}, {1778.8949, 1.0406552}, {1763.8195, 1.0392634}, {1748.7441, 1.0378224}, {1733.6688, 1.0363292}, {1718.5934, 1.0347802}, {1703.518, 1.0331707}, {1688.4426, 1.0314941}, {1673.3673, 1.0297415}, 
    {1658.2919, 1.0279002}, {1643.2165, 1.025952}, {1628.1412, 1.0238707}, {1613.0658, 1.0216182}, {1597.9904, 1.0191391}, {1582.915, 1.0163525}, {1567.8397, 1.0131417}, {1552.7643, 1.0093401}, {1537.6889, 1.0047152}, 
    {1522.6136, 0.99895508}, {1507.5382, 0.99166793}, {1492.4628, 0.9824059}, {1477.3874, 0.97072308}, {1462.3121, 0.95625418}, {1447.2367, 0.93878259}, {1432.1613, 0.91826926}, {1417.086, 0.89483803}, 
    {1402.0106, 0.86873404}, {1386.9352, 0.84027554}, {1371.8598, 0.80981198}, {1356.7845, 0.77769325}, {1341.7091, 0.74424993}, {1326.6337, 0.70978283}, {1311.5584, 0.67455941}, {1296.483, 0.63881506}, 
    {1281.4076, 0.60275746}, {1266.3322, 0.56657254}, {1251.2569, 0.53043096}, {1236.1815, 0.4944942}, {1221.1061, 0.45891987}, {1206.0307, 0.4238656}, {1190.9554, 0.3894916}, {1175.88, 0.35596148}, {1160.8046, 0.32344154}, 
    {1145.7293, 0.29209849}, {1130.6539, 0.26209576}, {1115.5785, 0.23358875}, {1100.5031, 0.20671934}, {1085.4278, 0.18161021}, {1070.3524, 0.1583594}, {1055.277, 0.1370358}, {1040.2017, 0.11767584}, 
    {1025.1263, 0.10028183}, {1010.0509, 0.084822106}, {994.97554, 0.07123293}, {979.90017, 0.05942197}, {964.8248, 0.049273035}, {949.74943, 0.040651617}, {934.67406, 0.033410795}, {919.59868, 0.027397081}, 
    {904.52331, 0.022455872}, {889.44794, 0.018436284}, {874.37257, 0.015195242}, {859.2972, 0.012600771}, {844.22182, 0.01053444}, {829.14645, 0.008892887}, {814.07108, 0.007588385}, {798.99571, 0.006548454}, 
    {783.92034, 0.005714659}, {768.84496, 0.005040865}, {753.76959, 0.004491226}, {738.69422, 0.004038206}, {723.61885, 0.003660784}, {708.54348, 0.003342951}, {693.46811, 0.003072499}, {678.39273, 0.002840079}, 
    {663.31736, 0.002638491}, {648.24199, 0.002462145}, {633.16662, 0.002306669}, {618.09125, 0.00216861}, {603.01587, 0.002045214}, {587.9405, 0.00193427}, {572.86513, 0.001833981}, {557.78976, 0.001742876}, 
    {542.71439, 0.00165974}, {527.63902, 0.001583562}, {512.56364, 0.001513494}, {497.48827, 0.001448819}, {482.4129, 0.001388928}, {467.33753, 0.001333298}, {452.26216, 0.00128148}, {437.18678, 0.001233084}, 
    {422.11141, 0.00118777}, {407.03604, 0.001145243}, {391.96067, 0.00110524}, {376.8853, 0.001067531}, {361.80992, 0.001031911}, {346.73455, 0.000998196}, {331.65918, 0.000966223}, {316.58381, 0.000935844}, 
    {301.50844, 0.000906924}, {286.43307, 0.00087934}, {271.35769, 0.00085298}, {256.28232, 0.00082774}, {241.20695, 0.000803522}, {226.13158, 0.000780232}, {211.05621, 0.000757781}, {195.98083, 0.00073608}, 
    {180.90546, 0.000715041}, {165.83009, 0.00069457}, {150.75472, 0.00067457}, {135.67935, 0.000654928}, {120.60397, 0.000635516}, {105.5286, 0.000616174}, {90.453231, 0.000596693}, {75.377859, 0.000576785}, 
    {60.302487, 0.000556011}, {45.227116, 0.000533637}, {30.151744, 0.000508209}, {15.076372, 0.000475883}, {0, 0.000410148}}; 

FreeElectronFraction := Interpolation[FreeElectronFractionData, InterpolationOrder -> 1]

ElectronNumberDensity[\[Eta]_] := (redShift = 6.64*^18^2/((c - Sqrt[c]*Sqrt[c - 2.*A*\[Eta]])/A)^2 - 1.; FreeElectronFraction[redShift]*1.42*^-7*(1. + redShift)^3)

Plot[NIntegrate[ElectronNumberDensity[eta], {eta, \[Eta], 3.78*^18}, MaxRecursion -> 15], {\[Eta], 1.47*^17, 2.66*^17}]

ListPlot[FreeElectronFractionData]

The 25 seconds doesn't seem to be a lot, but this calculation is inside another integral which didn't complete in eight hours. As nearly as I can tell, this integral is the culprit. Specifically the interpolation function.

I've seen other suggested solutions on this board, but none of them worked for me. One of the solutions looked promising: creating a pure function based on the interpolated data and using that in the integral, but that is beyond my skills.

Improve this question
$\endgroup$
2
  • $\begingroup$ You could try turning off symbolic processing and reduce the MaxRecursion NIntegrate[ElectronNumberDensity[eta], {eta, \[Eta], 3.78*^18}, Method -> {Automatic, "SymbolicProcessing" -> 0}, MaxRecursion -> 2]. $\endgroup$
    – Tim Laska
    Sep 4, 2019 at 1:44
  • $\begingroup$ Sorry, but reducing the MaxRecursion appears to introduce unacceptable inaccuracies to the integral results. The 'SymbolicProcessing' did seem to have a positive effect on performance, though. Thanks! $\endgroup$
    – Quark Soup
    Sep 4, 2019 at 12:15

2 Answers 2

Reset to default
14
$\begingroup$

The code in the question takes about 40 seconds on my computer. Turning off SymbolicProcessing, as suggested by Tim Laska in a comment, reduces the run time to about 16 seconds. Although reducing MaxRecursion further reduces the run time, it also generates error messages that the integration is not converging to sufficient accuracy. (Whether this matters from a practical point of view depends on how the results will be used.)

Instead, let us take advantage of the fact that NIntegrate in Plot repeatedly performs the same integral over the range {2.66*^17, 3.78*^18} to break the integral into two pieces.

int = NIntegrate[ElectronNumberDensity[eta], {eta, 2.66*^17, 3.78*^18}]
Plot[int + NIntegrate[ElectronNumberDensity[eta], {eta, η, 2.66*^17}, 
    Method -> {Automatic, "SymbolicProcessing" -> 0}], {η, 1.47*^17, 2.66*^17}, 
    ImageSize -> Large, LabelStyle -> {Bold, Black, 15}]

which produces the desired plot in under 5 seconds.

enter image description here

The plot agrees well with that produced by the code in the question.

Addendum

Even faster is to use NDSolveValue instead of NIntegrate.

NDSolveValue[{g'[eta] == -ElectronNumberDensity[eta], g[2.66*^17] == int}, 
    g, {eta, 1.47*^17, 2.66*^17}]
Plot[%[η], {η, 1.47*^17, 2.66*^17}, ImageSize -> Large, LabelStyle -> {Bold, Black, 15}]

which produces the same curve in less than a second.

Improve this answer
$\endgroup$
6
  • $\begingroup$ This looks promising, but I'm unfamiliar with the workings of ODE solving. Why, for instance, do you take $g'[\eta] == -ElectronNumberDensity[\eta]$ for the first equation? $\endgroup$
    – Quark Soup
    Sep 4, 2019 at 12:18
  • 4
    $\begingroup$ @Quarkly Basically, we are applying the fundamental theorem of calculus, that a differential is the inverse of an integral. Solving the differential equation, therefore, gives the integral, if boundary conditions are chosen properly. The minus sign arises, because η is the lower bound of the original integral, rather than the upper bound. Using NDSolve is much more efficient than using NIntegrate, when one endpoint of the integral is a variable. This can be seen by understanding that both functions in effect sum a set of finite elements, but only NDSolve returns a running total. $\endgroup$
    – bbgodfrey
    Sep 4, 2019 at 14:46
  • $\begingroup$ Can you tell me how you picked the lower and upper limits on the first NIntegrate function? Why did you pick the upper boundary of the Plot function call (2.66*^17) as the lower limit of the NIntegrate function? The boundary of the Plot is arbitrary. $\endgroup$
    – Quark Soup
    Sep 6, 2019 at 20:01
  • $\begingroup$ I chose the upper bound of the Plot as the lower bound for computing int, because doing so avoided recomputing the same quantity again and again. However, I did this before considering the NDSolve` approach. I could just as well set int == 0 and solved the ODE all the way to 3.78*^18, which is not noticeably slower.. $\endgroup$
    – bbgodfrey
    Sep 6, 2019 at 21:12
  • 1
    $\begingroup$ Use intg = NDSolveValue[{g'[etaprime] == -ElectronNumberDensity[etaprime], g[const2] == g2}, {etaprime,const1,const2}], where const1 is the smallest value of eta in which you are interested, and const2 is the largest value in which you are interested. Of course, you need to know g2, the value of the integral at const2. If the integrand becomes small beyond const2, you can set g2 == 0. intg then is the integral you want as an InterpolationFunction, and you can use intg[eta] however you wish. For instance, Plot[intg[eta], {eta, const1, const2}]. $\endgroup$
    – bbgodfrey
    Sep 7, 2019 at 1:49
8
$\begingroup$

This

c = 2.99792*10^5;
A = 3.87624*10^-14;
FreeElectronFractionData = {...YourData...}; 
Clear[\[Eta]]; (* just in case there was a prior assignment to Eta *)
redShift = 6.64*^18^2/((c - Sqrt[c]*Sqrt[c - 2.*A*\[Eta]])/A)^2 - 1.;
FreeElectronFraction = Interpolation[FreeElectronFractionData, InterpolationOrder -> 1];
ElectronNumberDensity[\[Eta]_] := FreeElectronFraction[redShift ]*1.42*^-7*(1. +redShift^3);
Plot[NIntegrate[ElectronNumberDensity[eta], {eta, \[Eta], 3.78*^18}, MaxRecursion -> 15], {\[Eta], 1.47*^17, 2.66*^17}]
ListPlot[FreeElectronFractionData]

completes in a less than ten seconds, mostly I believe because I changed

FreeElectronFraction := Interpolation[...]

to

FreeElectronFraction = Interpolation[...]

which means the calculation of the interpolating function will be done only once and the result will then be repeatedly used instead of calculating the interpolation function thousands of times.

I also moved the location of your calculation of redShift which may or may not have been necessary, but I was getting recursion errors and fixed those by moving the redShift calculation.

Please check all this very carefully to make certain the result is correct.

Improve this answer
$\endgroup$
6
  • 1
    $\begingroup$ Your \[Eta] might not being used properly. You should probably turn redShift into a function of \[Eta] and use that in ElectronNumberDensity. That way there's no question the order of operations is right. $\endgroup$
    – b3m2a1
    Sep 4, 2019 at 4:37
  • $\begingroup$ @b3m2a1 I was trying to do the minimum number of changes to his code, but you make an excellent point. If the original poster tests this carefully then I would hope that he would uncover errors. If my changes are incorrect then as suggested try redShift[\[Eta]_]:= 6.64*^18^2/((c-Sqrt[c]*Sqrt[c-2.*A*\[Eta]])/A)^2-1.; and turn the two uses of redShift into redShift[\[Eta]] and see if this corrects the problem. Thank you $\endgroup$
    – Bill
    Sep 4, 2019 at 5:05
  • $\begingroup$ The change in your comment indeed corrects the problem, but it also eliminates much of the run time savings from your answer. $\endgroup$
    – bbgodfrey
    Sep 4, 2019 at 5:07
  • $\begingroup$ Redshift is a function of conformal time (eta). Evaluating it once give the wrong answer. However, I'm surprised that this simple calculation causes such an impact on timing. $\endgroup$
    – Quark Soup
    Sep 4, 2019 at 12:16
  • 1
    $\begingroup$ @Quarkly the issue was that you used := for the interpolation. That causes the interpolation to be re-evaluated hundreds or thousands of time. Just create the interpolating function once. $\endgroup$
    – b3m2a1
    Sep 4, 2019 at 23:14

Your Answer

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

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