Cubic spline Interpolation

Question

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]

Answer 1

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

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.

Answer 2

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