2
$\begingroup$

I'd like to fit a cubic spline object to some data in mathematica. I think I'm ok in how to set it up, I just need a little help with some of the syntax. I'd like to generate a list of data first comprising about 100,000 points. The way I was doing this was to use a For loop:

For[i = 1, i <= 10, i++, data = {{i,  4*Pi/9/Log[i/0.09]}}]
{{10, 0.296413}}

But this only generated the element of data I wanted for i = 10 (for i up to 100000 the run time seems long, but to express my problem the number of points shouldn't matter). Could someone tell me how to fix the For loop so I can obtain 'data' as a brace containing the 10 sets of data points?

I'm also interested in creating a data set with points of the form e.g {{1.1, f(1.1)}, {1.2, f(1.2)}, {1.3, f(1.3)}, ... }, i.e., where the increment of my x values is 0.1 and not 1 as in the For Loop. Is there a way to modify the For Loop in this case too?

In the end I expect to evaluate something like

cubicsp1 = Interpolation[data]
$\endgroup$

closed as off-topic by Alexey Popkov, LCarvalho, MarcoB, m_goldberg, gwr Nov 24 '17 at 22:40

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question arises due to a simple mistake such as a trivial syntax error, incorrect capitalization, spelling mistake, or other typographical error and is unlikely to help any future visitors, or else it is easily found in the documentation." – Alexey Popkov, LCarvalho, MarcoB, m_goldberg, gwr
If this question can be reworded to fit the rules in the help center, please edit the question.

  • $\begingroup$ do not use For ( ever ). Look up Table .. data=Table[{i, 4*Pi/9/Log[i/0.09]},{i,10}] $\endgroup$ – george2079 Nov 21 '17 at 16:39
  • $\begingroup$ I did something along these lines here let me know if you need more. $\endgroup$ – Hugh Nov 21 '17 at 16:56
  • 1
    $\begingroup$ Note that Interpolation by default doesn't construct a smooth natural cubic spline if that's what you are going for, but some kind of piecewise (hermite?) spline (meaning that the first derivative is in general not continuous). $\endgroup$ – Thies Heidecke Nov 21 '17 at 18:31
  • $\begingroup$ @CAF The title doesn't reflect the essence of the question. Please make more adequate title. $\endgroup$ – Alexey Popkov Nov 22 '17 at 8:09
2
$\begingroup$

There are much easier ways to construct the data than a For loop. george2079 shows one in the comments, here is another that uses Map (represented by /@) and Slot (represented by #), which essentially fills all the slots in the specified Range. Then f is defined to be an interpolating function (order 3 specifies cubic), and you can evaluate the function f at any point in the range.

dat = {#, 4*Pi/9/Log[#/0.09]} & /@ Range[1, 10, 0.1];
f = Interpolation[dat, InterpolationOrder -> 3];
{f[2.2], f[2.211], f[5.03]}
$\endgroup$
  • $\begingroup$ Thanks! I should say I'm also interested in using the interpolating function I get from this analysis as an integrand involved in a numeric integral. Calling NIntegrate[cubicsp1[j],{j,1,100000}] seems to be taking a long time. Could you perhaps edit your answer (if you know a way) to modify this line or the way I've set up the interpolating method so that the numerical integration is done quicker? $\endgroup$ – CAF Nov 24 '17 at 18:23
  • $\begingroup$ You started out with data, then created an interpolating function, now you want to integrate? Why not just sum the data in the first place? The interpolating function does not provide any new accuracy over what is in the data itself. $\endgroup$ – bill s Nov 24 '17 at 20:10
  • $\begingroup$ Ah sorry in order to express the problem I was having with the interpolation, I seemed to have oversimplified what I'm trying to do - basically I have an integral of the form $\int_1^{15000} f_1(k) f_2(k) dk$ where $f_1(k)$ is the 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}]; and then f1 = Interpolation[data];, where RunDec is a specialised mathematica package - I didn't mention this in OP because of its rarity onsite $\endgroup$ – CAF Nov 24 '17 at 21:10
  • $\begingroup$ Maybe the problem lies with your (unspecified) f2. If I run: dat = {#, 4*Pi/9/Log[#/0.09]} & /@ Range[1, 100000, 0.1]; f = Interpolation[dat, InterpolationOrder -> 3]; NIntegrate[f[x], {x, 1, 100000}] // Timing it returns in 1.3 seconds. $\endgroup$ – bill s Nov 25 '17 at 0:59
2
$\begingroup$

It wasn't quite clear if you wanted 100,000 points between 1 and 10 or 100,000 points from one to 100,000.

I will assume the former (changing it to the latter is trivial).

I like to use Table rather than Map because the syntax is easier to follow and the performance almost as good in most cases.

The iterator {i, 1, 10, 1/10000} creates 100,000 points between 1 and 10.

data = Table[{i, (4 π)/(9 Log[(10 i)/9])}, {i, 1, 10, 1/10000}];

Create the interpolation function

f = Interpolation[data, InterpolationOrder -> 3];

Compare the data and function

Show[
 ListPlot[data, PlotRange -> All, PlotStyle -> Black],
 Plot[f[x],  {x, 1, 10}, PlotStyle -> {Red, Dashed}]
 ]

Mathematica graphics

On this scale they match perfectly.

Part 2 - NIntegrate

I don't know how many times you are doing this but I found NIntegrate to be reasonable with the interpolating function.

NIntegrate[f[x], {x, 1, 10}] // AbsoluteTiming
(* {0.546858, 10.3072} *)

It took about 0.5 seconds on my system (Mathematica 11.1.1, Windows 7, 64 bit, 3.5GHz, 16.0 GB RAM).

In order to speed it up you might reduce the granularity of the data input to the interpolating function. Also, since it becomes very smooth and flat beyond 3 or 4 seconds, make a further reduction of the density in that region.

For for example, below the density is reduced by a factor of two from 1 to 3, and by a factor of 10, from 3 to 10.

dataReduced = Join[
   Table[{i, (4 π)/(9 Log[(10 i)/9])}, {i, 1, 3, 2/10000}],
   Table[{i, (4 π)/(9 Log[(10 i)/9])}, {i, 3 + 1/1000, 10, 1/1000}]
   ];

 fR = Interpolation[dataReduces, InterpolationOrder -> 3];

NIntegrate[f[x], {x, 1, 10}] // RepeatedTiming
(* {0.53, 10.3072} *)

NIntegrate[fR[x], {x, 1, 10}] // RepeatedTiming
(* {0.12, 10.3072} *)

This results in a speedup of about a factor of 4 - 5.

Another approach would be to make an Interpolation of the integral.

Your input function has a closed form for the integral.

Integrate[4 π/(9 Log[(10 x)/9]), x]
(* 2/5 π LogIntegral[(10 x)/9] *)

Use that to make a Table and an Interpolation.

dataIntegral = Table[{i, 2/5 π LogIntegral[(10 i)/9]}, {i, 1, 10, 1/10000}];

fI = Interpolation[dataIntegral];

and then

fI[10] - fI[1] // RepeatedTiming // N
(* {0.0000709202, 10.3072} *)

This is blindingly fast.

The answer is particular to your input function which was a closed form for the integration.

What if your input function does not have a closed form? Say it is an arbitrary function fInput.

Then you could pay the price once as follows:

dataIntegral = Table[{i, NIntegrate[fInput[x], {x, 1, i}]}, {i, 1, 10, 1/10000}];

fI = Interpolation[dataIntegral]

That might take a long time to create but once done you will have a rapid response when you invoke fI.

$\endgroup$
  • $\begingroup$ Thanks! I should say I'm also interested in using the interpolating function I get from this analysis as an integrand involved in a numeric integral. Calling NIntegrate[cubicsp1[j],{j,1,100000}] seems to be taking a long time. Could you perhaps edit your answer (if you know a way) to modify this line or the way I've set up the interpolating method so that the numerical integration is done quicker? $\endgroup$ – CAF Nov 24 '17 at 17:19
  • $\begingroup$ I'm new to this way of doing things but I can guess it's something to do with the integral looping over the table of values at every value it needs which might explain the longevity of the running time. $\endgroup$ – CAF Nov 24 '17 at 17:25
  • $\begingroup$ No problem! I have done so, it is here mathematica.stackexchange.com/questions/160622/… if you are able to help out. Thanks $\endgroup$ – CAF Nov 24 '17 at 21:54

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