Interpolating Functions

Introduction

This notebook explains how to use and operate on interpolating functions in Mathematica. These functions are either created by certain commands, like Interpolation, or else are returned by NDSolve as the solutions to differential equations. In both cases the functions are intended to predict continuous behavior from a discrete set of data points, and to behave as continuous functions with respect to most other Mathematica commands. There are, however, subtle differences between NDSolve interpolations and Interpolate interpolations, which are (this is Mathematica, after all) very frustrating to figure out. This notebook, then, is also intended to explain these differences.

Plain Old Interpolation

The simplest way to creat interpolating functions in Mathematica is to use the command Interpolation. This command accepts as an argument an array of data points. The array may consist of either values or argument-value pairs (i.e., either {y1, y2, y3...} or {{x1, y1},{x2, y2}...}). Interpolation examines the array and predicts continuous behavior from it - in other words, it interpolates all values falling between those given in the array. It is very powerful and very robust. Here is an example.

Creating Interpolating Functions

First, make a table of sample data. For simplicity, use a sine function:

data = Table[{t,Sin[ t]},{t,0, 4Pi+.1, .1}];
ListPlot[data]

Now apply Interpolation to the data array:

dataIF = Interpolation[data]
Plot[dataIF[t],{t,0,4Pi}]

Notice how, once dataIF is set equal to the output of Interpolation, dataIF behaves as a function - you can add a [t] to the end to evaluate it at specific points, or plot it. 
In fact, you can do almost anything to dataIF that you would do to an ordinary function. Examples:

FindMaximum[dataIF[t],{t,Pi/2}]

FindRoot[dataIF[t],{t,Pi}]

However, if you wish to set the output of functions like FindMaximum or FindRoot equal to some variable, and then operate on that variable, you must be very careful. The output of FindMaximum or FindRoot is technically an array of rules. To get the actual value, you must apply the rule contained in the element of the array to some variable. Here's an example.

root = FindRoot[dataIF[t],{t,Pi}]

Operating on root won't work:

value = dataIF[root]

You must take the first (and only) element of the array that is actually returned by FindRoot, which is a rule, apply the rule to the given variable, and set the result equal to a different variable  (or the same one, if you don't mind overwriting it). Rules are essentially just Mathematica's way of specifying transformations. They always use the arrow notation, and can be applied to a variable or expression using the /. combination. Read /. as "such that" and read → as "goes to". For example:

arule = a -> (x+1);
a^2 /. arule

can be read as " (a^2) such that a goes to x+1".

However, it is important to note that the variable used in the declaration of the rule must be present in the expression the rule is applied to. For example, the following code won't work:

arule = a -> (x+1);
b^2 /. arule

Extracting the output of FindRoot is very similar to the above procedure. Use the First command to get the element of the array, and apply this element (which is a rule) to some variable:

root

First[root]

troot = t /. First[root];
value = dataIF[troot]

troot

Note how now dataIF, or any function, will accept the variable troot as an argument, since it is no longer an array but an actual number.

Taking Derivatives

Taking derivatives of interpolating functions is quite straightforward. It requires using the Derivative command, which is more general than the usual command for symbolic expressions, D[]. The output is simply another interpolating function. Here's an example:

dataIFd = Derivative[1][dataIF]
Plot[{dataIF[t], dataIFd[t]},{t,0,12}]

Note that you do not need to specify an argument (i.e., no [t]) when applying Derivative to dataIF. The reasons are obscure and confusing, having to do with how function objects work in Mathematica. The [1] in brackets after Derivative just specifies how many times to take the derivative - 2 or 3 will give the 2nd or 3rd derivative, and so on.

Integrals

Unfortunately, you can't use the Integrate command to find an interpolating function approximating the symbolic integral of a given interpolating function, as you can with the Derivative command. However, you can find numerical integrals very easily. Here's an example:

definiteInt = Integrate[dataIF[t],{t,0,Pi}]

If you don't mind a little bit of extra code, it's possible to manually create an interpolating function that approximates the symbolic integral, using the Table command. The resulting function isn't centered on the x-axis, which is a consequence of the arbitrary integration constant produced during indefinite integration, but the 'shape' of the function is pretty well approximated.

dataIFi = Interpolation[Table[{i,Integrate[dataIF[t],{t,0,i}]},{i,0,12,.1}]]

Plot[dataIFi[t],{t,0,12}]

Notes on Power

The Interpolation command is very powerful. This is illustrated below, with the same sine wave as used above, except with drastically reduced sampling. The first plot shows the sampled data points, of which there are only six. The second plot, of the interpolation function produced from these points, superimposed on the original sine wave, illustrates how close to the original function the interpolation is.

data2 = Table[{t,Sin[ t]},{t,0, 4Pi+.1, 2}];
ListPlot[data2,PlotStyle->{PointSize[.015]}]
Plot[{Sin[t],Interpolation[data2][t]},{t,0,12}]

NDSolve Interpolating Functions

The other important use of interpolating functions in Mathematica is as the output of the NDSolve command. Essentially, NDSolve generates a list of solutions to the given differential equations in the form {t, x(t)}, where x(t) is the solution for some dynamical variable. NDSolve then automatically applies Mathematica's Interpolation command to this list, and outputs the resulting interpolating function. Here's a simple example.

Creating Interpolating Functions

sol = NDSolve[{x''[t] - .2x'[t] + 5x[t] ==0, y'[t] == x[t]+.5,x[0] ==1, z'[t] == .5,x'[0] == 10 ,y[0] ==0,z[0] ==0},{x,y,z},{t,0,100}]

There are several important things to note about the form of NDSolve's output. First, it is a nested array - there is the top-level array, the element of which is a two element array containing the interpolating functions. Second, the interpolating functions are given in terms of rules, so that you must use the /. notation to 'access' the interpolating functions. First gives the first element in the expression. For example the first element in this list is:

First[{a,b,c}]

while the first element in this list is a list :

First[{{a,b},{d,c}}]

First[{{a,b,c}}]

and the only element in this list is a list containing three elements :

First[sol]

First[sol] is a list of three rules

First[First[sol]]

First[First[sol]] is the first rule, or equivalently

Part[First[sol],1]

Part[First[sol],2]

y[t] /.sol

y[t] /. sol is a list of one element, the function y[t]

y[t] /.First[sol]

y[t] /. First[sol] is the function y[t]

y[10] /.First[sol]

y[10] /.sol

First[y[10] /.sol]

First [y[t] /. sol] is also the function y[t]

The command FindRoot can be used exactly as they are with regular interpolating functions:

FindRoot[x[t] ==0 /. sol,{t,1}]

Plot[x[t]/.sol,{t,0,5}]

FindMaximum works too:

FindMaximum[x[t]/.sol,{t,3}]

Creating Plots

Plotting the interpolating functions produced by NDSolve is not as straightforward as plotting regular interpolating functions. The difficulties again arise from one-element arrays and quirks in the plotting functions. 
Plotting one variable is, however, easy:

Plot[x[t]/.sol,{t,0,20}]

To plot two variables on the same graph:

Plot[{x[t],y[t]} /. sol,{t,0, 20}]

ParametricPlot[{x[t],y[t]} /. sol,{t,0,20}]

ParametricPlot3D[{x[t],y[t],z[t]} /. sol, {t,0,20}]

Getting actual data from Interpolating Functions

There are some commands in Mathematica that need lists or arrays of data, and not interpolating functions. An important example are the commands for doing Fourier transforms. While Mathematica does have commands for performing Fourier transforms on analytic expressions, they do not work with interpolating functions. It is therefore necessary to create a list or array of data from an interpolating function, and use the discrete Fourier transform command on this data. This is straightforward to do, using the Table command. The only trick is to once again use the First command. Here is an example, using the interpolating functions used above.

data = Table[First[x[i] /. sol],{i,0,5,.1}]

data = Table[First[x[i] /. sol],{i,0,50,.02}];
dataF = Fourier[data];
ListPlot[data,PlotRange->All]
ListPlot[{Abs[dataF],Abs[dataF]},Joined->{True,False},PlotRange->{{0,100},All}]