Interpolation

Given a set of $x$ and $y$ points, estimate intermediate values.
This is normally done be fitting a function to local points and then using the function to evaluate the desired points.
Example, given a set of points on a nonuniform grid, interpolate values to a uniform grid.
Interpolation will deal with 4 sets of data:
- Given: $x_{\text{given}}$, $y_{\text{given}}$, $x_{\text{desired}}$
- Find: $y_{\text{desired}}$.

Linear interpolation

For a given $x_d$ (“d” for “desired”), you need to find your neighborhood.
- Find the two nearest values of the given data $x_g$ that bracket $x_d$.
- Really, you are finding the two list indicies of $x_g$ that bracket $x_d$.
- Call these indices $i_{lo}$ and $i_{hi}$.
The equation for the line through the given bracketing points is

$$y = y_g[i_{lo}] + (x - x_g[i_{lo}])\frac{y_g[i_{hi}]-y_g[i_{lo}]}{x_g[i_{hi}]-x_g[i_{lo}]}.$$

Evaluate this at $x=x_d$ to get $y_d$:

$$y_d = y_g[i_{lo}] + (x_d - x_g[i_{lo}])\frac{y_g[i_{hi}]-y_g[i_{lo}]}{x_g[i_{hi}]-x_g[i_{lo}]}.$$

Extrapolation

If $x_d<x[0]$, then $i_{lo}=0$ and $i_{hi}=1$.
If $x_d>x[n-1]$, (for $n$ given points), then $i_{lo}=n-2$ and $i_{hi}=n-1$.
The above formula still works, but $x_d$ is not bracketed. We call this extrapolation.

Python

Python has a built in interpolator called interp1d.
This will interpolate a 1-D array.

Example

Take given x data as 10 linearly spaced points between 0 and 10.
Take given y data from the function $\cos(x^2/8)$.

import numpy as np
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
%matplotlib inline

def plot_a():
    x_g= np.linspace(0,10,10)         # Interpolate between the given data
    y_g= np.cos(x_g**2.0/8.0)      
    xx      = np.linspace(0,10,1000)  # for plotting
    yy      = np.cos(xx**2.0/8.0)     # for plotting
    
    x_i     = np.linspace(0,10,100)   # interpolate to these points.
    
    #------- Linear interpolation  
    
    f_linear   = interp1d(x_g, y_g)   # get an interpolation function
    y_i_linear = f_linear(x_i)        # interpolate to y_i at x_i
    
    #------- Plot the results
    
    plt.plot(x_g, y_g,'o')
    plt.plot(x_i, y_i_linear, '-')
    plt.plot(xx, yy, ':')
    plt.legend(['data', 'linear', 'perfect'], loc='best', frameon=False)
    plt.ylim([-1.5,1.5])
plot_a()

Cubic Splines

Cubic splines give a higher order local fit to the data points.
Fit a cubic polynomial for region $i$ to two adjacent data points $x_i$ and $x_{i+1}$ $$p_i(x) = a_i + b_i x + c_ix^2 + d_ix^3.$$
Suppose we have a list of 5 points $x_0$ through $x_4$, with corresponding $y_0$ through $y_4$.
There are four intervals on which splines are defined.

def plot_b():
    x_g = np.linspace(0,4,5)
    y_g = np.array([1,3,2.5,0.5,1])
    fs  = interp1d(x_g,y_g, kind='cubic')
    
    plt.rc('font', size=12);
    fig, ax = plt.subplots(1,1);
    for i in range(4):
        xx = np.linspace(i,i+1,1000)
        ax.plot(xx, fs(xx), '-')
    ax.plot(x_g,y_g,'ko')
    ax.set_ylim([-0.1,3.5])
    
    for i in range(4):
        ax.text(i+0.4,0.5*(fs(i)+fs(i+1))-0.2,r'$I_'+str(i)+'$')
    for i in range(5):
        ax.text(i-0.1,0.2,r'$x_'+str(i)+'$')
    for i in range(5):
        ax.text(i-0.1,fs(i)+0.25,r'$y_'+str(i)+'$')
plot_b()

Again, we have 5 points and 4 splines.
Each spline has 4 coefficients, so we have 16 degrees of freedom.
Need 16 constraints to specify the splines and solve for the spline polynomial coefficients.
- Each spline passes through two points: 8
- Let the derivative $p^{\prime}(x)$ be continuous at the spline interfaces: 3
- Let the second derivative also be continuous at the spline interfaces: 3
- This gives 14 constraints. We need two more.
  - Specify $p^{\prime\prime}(x)$ at the two ends: 2.
This results in a linear system of equations that can be solved directly.
Hoffman and Press et al. (Numerical Recipies) reformulate the equations in a manner that gives a tridiagonal system of equations.
Note that the whole domain is coupled when solving the splines, and all of the splines are solved together.

Equations

Hoffman p. 221 and Numerical Recipes p. 120 give details of the equations.

A general cubic equation is given by $$f(x) = ax^3 + bx^2 + cx + d.$$ The second derivative is $$f^{\prime\prime}(x) = 6ax + 2b,$$ which is a linear profile. Hence, the second derivative profile can be written as $$f^{\prime\prime} = Af_i^{\prime\prime} + Bf_{i+1}^{\prime\prime},$$ where $$A = \frac{x_{i+1}-x}{x_{i+1}-x_i},$$ $$B = 1-A = \frac{x-x_i}{x_{i+1}-x_i}.$$

If we integrate this twice, we get an expression for $f(x)$ and $f^\prime(x)$, with two constants of integration. These constants are evaluated by requiring that $f(x) = f(x_i)$ at the two end points of a given spline. This results in the following spline equation:

$$f(x) = \frac{f_i^{\prime\prime}}{6(x_{i+1}-x_i)}(x_{i+1}-x)^3 + \frac{f_{i+1}^{\prime\prime}}{6(x_{i+1}-x_i)}(x-x_i)^3 + \\ \left[\frac{f_i}{x_{i+1}-x_i} -\frac{f_i^{\prime\prime}(x_{i+1}-x_i)}{6}\right](x_{i+1}-x) + \\ \left[\frac{f_{i+1}}{x_{i+1}-x_i} -\frac{f_{i+1}^{\prime\prime}(x_{i+1}-x_i)}{6}\right](x-x_i)$$

This equation is used to evaluate the spline. But it is in terms of the second derivatives at the node points, which need to be evaluated. That evaluation is done by equating the first derivative between neighboring splines. This yeilds a tridiagonal system for the second derivatives:

$$(x_i - x_{i-1})f_{i-1}^{\prime\prime} + 2(x_{i+1}-x_{i-1})f_i^{\prime\prime} + (x_{i+1}-x_i)f_{i+1}^{\prime\prime} = 6\frac{f_{i+1}-f_i}{x_{i+1}-x_i} - 6\frac{f_i-f_{i-1}}{x_i-x_{i-1}},$$

which can be solved using the Thomas Algorithm.

Note, the boundary values of $f^{\prime\prime}$ are needed, and these are usually set to zero.

Python

For cubic splines, just add the keyword: kind=‘cubic’ to the interp1d function call.

def plot_c():
    x_g= np.linspace(0,10,10)         # Interpolate between the given data
    y_g= np.cos(x_g**2.0/8.0)      
    xx      = np.linspace(0,10,1000)  # for plotting
    yy      = np.cos(xx**2.0/8.0)     # for plotting
    
    x_i     = np.linspace(0,10,100)   # interpolate to these points.
    
    #------- Linear interpolation  
    
    f_cubic   = interp1d(x_g, y_g, kind='cubic')   # get an interpolation function
    y_i_cubic = f_cubic(x_i)        # interpolate to y_i at x_i
    
    #------- Plot the results
    
    plt.plot(x_g, y_g,'o')
    plt.plot(x_i, y_i_cubic, '-')
    plt.plot(xx, yy, ':')
    plt.legend(['data', 'cubic spline', 'perfect'], loc='best', frameon=False)
    plt.ylim([-1.5,1.5])
plot_c()

Find your location

To apply linear interpolation at a given point, we first need to find the location of $x_d$ in the list.
- That is, we need $i_{lo}$ and $i_{hi}$.
Interpolation routines will assume the $x$ data are monotonically increasing or decreasing, and will issue an error otherwise.
Two cases:
1. Uniformly spaced points,
  - (or some other analytic spacing)
2. Arbitrarily spaced points.

Case 1: Uniform points

def locate_uniform(x, xd):
    Δx = x[1]-x[0]
    ilo = int((xd-x[0])/Δx)
    ihi = ilo+1
    if ilo < 0:
        return 0,1
    if ihi >= len(x):
        return len(x)-2, len(x)-1
    return ilo, ihi

x = np.linspace(2,5, 7)

for i in range(len(x)):
    print(f"{i}, {x[i]:.1f}")

xd = 3.7
ilo, ihi = locate_uniform(x, xd)

print(f"\nilo={ilo}, ihi={ihi} bracket point x={xd}")

0, 2.0
1, 2.5
2, 3.0
3, 3.5
4, 4.0
5, 4.5
6, 5.0

ilo=3, ihi=4 bracket point x=3.7

Case 2: Arbitrarily spaced points

Many approaches are possible.
Press et al. advocate using bisection.
Finds location in $\log_2N$ tries.

x = np.random.rand(20)
x = np.sort(x)

for xx in x:
    plt.plot([xx,xx],[0.4,0.6], '-', lw=0.5, color='gray')
    plt.plot(xx,0.5, 'o', ms=2, color='black')
plt.ylim([0,1])
plt.yticks([])
plt.xticks([]);

array([0.06307574, 0.07106583, 0.15084623, 0.19699816, 0.22199305,
       0.22368445, 0.36857857, 0.38526137, 0.38674107, 0.43815616,
       0.54737735, 0.5603345 , 0.61735068, 0.6294096 , 0.69707346,
       0.77795222, 0.8417913 , 0.85636799, 0.94958611, 0.97282627])

def getLoc(xgrid, x):                # bisection, from Numerical Recipes
    n = len(xgrid)
    jl, ju = 0, n-1
    while ju-jl > 1:
        jm = (ju + jl) >> 1;      # same as int((ju+jl)/2)
        if x >=xgrid[jm]:
            jl = jm
        else:
            ju = jm
    return jl, ju

getLoc(x, 0.5)

$\displaystyle \left( 9, \ 10\right)$

Example

Search the following grid for the bounding indices of $x=0.75$. The algorithm has the following progression, shown schematically below.

$j_l$	$j_u$	$j_m$	Condition	new point
0	6	3	$0.75\ge x_g[3]=0.4$	$j_l=3$
3	6	4	$0.75\ge x_g[4]=0.5$	$j_l=4$
4	6	5	$0.75\ge x_g[5]=1.0$	$j_u=5$
4	5	–	$j_u-j_l=1\rightarrow$ done

x = np.array([0.1,0.2,0.3,0.4,0.5,      1.0, 2.0])

fig, ax = plt.subplots(1,1, figsize=(10,3)) 

for xx in x:
    ax.plot([xx,xx],[0.6,0.8], '-', lw=0.5, color='gray')
    ax.plot(xx,0.7, 'o', ms=2, color='black')
    ax.plot([0.1, 2.0], [0.4,0.4], 'o', ms=6, color='blue')
    ax.plot([0.4, 2.0], [0.3,0.3], 'o', ms=6, color='green')
    ax.plot([0.5, 2.0], [0.2,0.2], 'o', ms=6, color='orange')
    ax.plot([0.5, 1.0], [0.1,0.1], 'o', ms=6, color='red')
    #ax.plot([0.75,0.75],[0.5,0.9], '--', lw=1.5, color='gray')
    ax.plot([0.75],[0.7], '*', ms=8, color='gray')
ax.set_ylim([0,1])
ax.set_yticks([])
ax.set_xticks(x);