1. Interpolation

Interpolation is a technique to infer information from a discrete and limited set of data. For instance, it is the root of numerical integration schemes, where one infer the shape of the curve between discrete points following some interpolation pattern.

Here we’ll discuss simple interpolation schemes: linear and quadratic interpolation; and an introduction to quadratic spline interpolation. For more information, please check the book by Tao Pang and the Wikipedia: Interpolation.

As a simple example, we’ll consider this figure of sine function with only 5 points in all cases below.

../_images/sin-5pts.png

1.1. Linear interpolation

A linear interpolation considers a straight line connecting each pair of consecutive points. The set of 5 discrete points above are \((x_i, y_i)\), with i=0..4. To interpolate we must fine the equation of aline connecting the dots. This is very simple, but let’s do it in details, since the same approach can be used for more complex interpolations.

We want to find the coefficients from \(y(x) = a x + b\) within the domain \([x_i, x_{i+1}]\). We use the known points to write the equations

\[\begin{split}\begin{align} y_i &= a x_i + b \\ y_{i+1} &= a x_{i+1} + b \end{align}\end{split}\]

These equations can be casted in a matrix form

\[\begin{split}\begin{pmatrix} x_i & 1 \\ x_{i+1} & 1 \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix} = \begin{pmatrix} y_i \\ y_{i+1} \end{pmatrix}\end{split}\]

This is a very simple equation that can be solved by hand, but for more complex cases you could used np.linalg.solve(…) to solve the linear system of equations. You’ll find

\[\begin{split}\begin{align} a &= \dfrac{y_{i+1} - y_i}{x_{i+1} - x_i} \\ b &= \dfrac{y_{i} x_{i+1} - y_{i+1} x_i}{x_{i+1} - x_i} \end{align}\end{split}\]

Applying this linear interpolation to the sine function above, we get this figure below.

Linear interpolation
The code
../_images/sin-linear.png 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 25})

#----------------------------
# my linear interp function
#----------------------------
def linear_interp(x, y, npts):
    # init arrays
    xs = np.linspace(x[0], x[-1], npts)
    ys = np.zeros_like(xs)
    # init first point
    ys[0] = y[0]

    # loop over sections set by original points
    for i0 in range(len(x)-1):
        # init line equation at this section
        a = (y[i0+1]-y[i0])/(x[i0+1]-x[i0])
        b = (y[i0]*x[i0+1]-y[i0+1]*x[i0])/(x[i0+1]-x[i0])
        # extract xs indexes within x-section range
        # i.e., points where x[i0] < xs < x[i0+1]
        js = np.argwhere((xs > x[i0]) & (xs <= x[i0+1]))[:,0]
        # apply linear interpolation over this range
        ys [js] = a*xs[js] + b
    # return results
    return xs, ys

# generate original points
x = np.linspace(0, 2*np.pi, 5)
y = np.sin(x)# + np.random.normal(0, 0.2, 7)
# exact sine with a lot of points as a reference
xe = np.linspace(0, 2*np.pi, 100)
ye = np.sin(xe)
# my interpolation over the 5 points above
x1, y1 = linear_interp(x, y, 25)

# plot all
plt.plot(xe, ye, ls=':', c='lightgray')
plt.scatter(x, y, c='black', s=100)
plt.scatter(x1, y1, c='red')
plt.xlabel(R'$\theta$ [rad]')
plt.ylabel(R'$\sin\theta$')
plt.xticks([0, np.pi/2, np.pi, 3*np.pi/2, 2*np.pi], ["0", R"$\dfrac{\pi}{2}$", R"$\pi$", R"$\dfrac{3\pi}{2}$", R"$2\pi$"])
plt.grid()
plt.tight_layout()
plt.show()

1.2. Quadratic interpolation

TO DO

1.3. Quadratic spline

TO DO

For more details check:

1.4. Numpy and scipy methods

TO DO