Math 441

Interpolation: Cubic Splines

What is a Spline?

Given \(n+1\) data points \((x_i,y_i)\) such that \(x_{i-1} < x_{i}\) for \(i = 1, 2, \dots, n\), spline of degree \(n\) is a function \(f:[x_0,x_n] \to \mathbb{R}\) such that

  • \(f(x_i) = y_i\) for \(i = 0, 1, \dots, n\),
  • \(f\) restricted to \([x_{i-1},x_i)\) is a polynomial of degree \(n\), for \(i = 1, 2, \dots, n\),
  • \(f\) satisfy some additional smoothness conditions.
  • Examples:
    • 0 degree splines: piecewise constant functions
    • linear splines: piecewise linear function, continuous at \(x_i\).
    • quadratic splines: piecewise quadratic function, continuously differentiable at \(x_i\).
    • cubic splines: piecewise cubic function, two times continuously differentiable at \(x_i\).

Requirements for Cubic Spline

  • \(S(x_i) = y_i\):
  • \(S\) is continuous
  • \(S'\) is continuous
  • \(S''\) is continuous

\(S_i(x_{i-1}) = y_{i-1}\)

\(S_i(x_{i}) = y_{i}\)

\(S'_i(x_{i}) = S'_{i+1}(x_i)\)

\(S''_i(x_{i}) = S''_{i+1}(x_i)\)

for \(i = 1, 2, \dots, n\)

for \(i = 1, 2, \dots, n\)

for \(i = 1, 2, \dots, n-1\)

for \(i = 1, 2, \dots, n-1\)


\[S(x) = \begin{cases} S_1 (x) & \text{if $x_0 \le x \lt x_1$}\\ S_2 (x) & \text{if $x_1 \le x \lt x_2$}\\ \quad\quad\vdots & \quad\quad\vdots \\ S_{n} (x) & \text{if $x_{n-1} \le x \le x_n$} \end{cases} \]

Coefficients

Write \(S_i\) as \[S_i (x) = {\color{green}{a_i}} + {\color{green}{b_i}}{\color{blue}{(x-x_{i-1})}} + {\color{green}{c_i}}{\color{blue}{(x-x_{i-1})}}^2 + {\color{green}{d_i}}{\color{blue}{(x-x_{i-1})}}^3\] Then

\[\begin{align*} S'_i (x) &= {\color{green}{b_i}} + 2{\color{green}{c_i}}{\color{blue}{(x-x_{i-1})}} + 3{\color{green}{d_i}}{\color{blue}{(x-x_{i-1})}}^2\\ S''_i (x) &= 2{\color{green}{c_i}} + 6{\color{green}{d_i}}{\color{blue}{(x-x_{i-1})}} \end{align*}\]

Define \(h_i = x_{i} - x_{i-1}\) for \(i = 1,\dots,n\). Then \[\begin{align*} S_i(x_{i-1}) &= {\color{green}{a_i}} & S_i(x_{i}) &= {\color{green}{a_i}} + {\color{green}{b_i}} {\color{blue}{h_i}} + {\color{green}{c_i}}{\color{blue}{h_i}}^2 + {\color{green}{d_i}}{\color{blue}{h_i}}^3 \\ S'_i(x_{i-1}) &= {\color{green}{b_i}} & S'_i(x_{i}) &= {\color{green}{b_i}} + 2{\color{green}{c_i}} {\color{blue}{h_i}} + 3{\color{green}{d_i}}{\color{blue}{h_i}}^2\\ S''_i(x_{i-1}) &= 2{\color{green}{c_i}} & S''_i(x_{i}) &= 2{\color{green}{c_i}} + 6{\color{green}{d_i}}{\color{blue}{h_i}} \end{align*}\]

Equations

\[\begin{align*} S_i(x_{i-1}) &= {\color{green}{a_i}} \\ S_i(x_{i}) &= {\color{green}{a_i}} + {\color{green}{b_i}} {\color{blue}{h_i}} + {\color{green}{c_i}}{\color{blue}{h_i}}^2 + {\color{green}{d_i}}{\color{blue}{h_i}}^3 \\ S'_i(x_{i-1}) &= {\color{green}{b_i}} \\ S'_i(x_{i}) &= {\color{green}{b_i}} + 2{\color{green}{c_i}} {\color{blue}{h_i}} + 3{\color{green}{d_i}}{\color{blue}{h_i}}^2\\ S''_i(x_{i-1}) &= 2{\color{green}{c_i}} \\ S''_i(x_{i}) &= 2{\color{green}{c_i}} + 6{\color{green}{d_i}}{\color{blue}{h_i}} \end{align*}\]

\[\begin{align*} S_i(x_{i-1}) &= y_{i-1} & \text{ \ \ for } i &= 1, 2, \dots, n\\ S_i(x_{i}) &= y_{i} & \text{ \ \ for } i &= 1, 2, \dots, n\\ S'_i(x_{i}) &= S'_{i+1}(x_i) & \text{ \ \ for } i &= 1, 2, \dots, n-1\\ S''_i(x_{i}) &= S''_{i+1}(x_i) & \text{ \ \ for } i &= 1, 2, \dots, n-1 \end{align*}\]


  • \(a_i = y_{i-1}\)
  • \(a_i + b_ih_i + c_ih_i^2 + d_ih_i^3 = y_{i}\)
  • \(b_i + 2c_ih_i + 3d_ih_i^2 = b_{i+1}\)
  • \(2c_i + 6d_ih_i = 2c_{i+1}\)
  • for \(i = 1, 2, \dots, n\)
  • for \(i = 1, 2, \dots, n\)
  • for \(i = 1, 2, \dots, n-1\)
  • for \(i = 1, 2, \dots, n-1\)

Simplifying

  • \(y_{i-1} + b_ih_i + c_ih_i^2 + d_ih_i^3 = y_{i}\) for \(i = 1, 2, \dots, n\)
  • \(b_i + 2c_ih_i + 3d_ih_i^2 = b_{i+1}\) for \(i = 1, 2, \dots, n-1\)
  • \(2c_i + 6d_ih_i = 2c_{i+1}\) for \(i = 1, 2, \dots, n-1\)

  • \(b_i + c_ih_i + d_ih_i^2 = \frac{y_{i} - y_{i-1}}{h_i}\) for \(i = 1, 2, \dots, n\)
  • \(b_i + 2c_ih_i + 3d_ih_i^2 = b_{i+1}\) for \(i = 1, 2, \dots, n-1\)
  • \(2c_i + 6d_ih_i = 2c_{i+1}\) for \(i = 1, 2, \dots, n-1\)

  • \(b_i = s_i - c_ih_i - d_ih_i^2\) for \(i = 1, 2, \dots, n\), where \(s_i = \frac{y_{i} - y_{i-1}}{h_i}\)
  • \(b_i + 2c_ih_i + 3d_ih_i^2 = b_{i+1}\) for \(i = 1, 2, \dots, n-1\)
  • \(d_i = \frac{2c_{i+1} - 2c_i}{6h_i}\) for \(i = 1, 2, \dots, n-1\)

Small change of variable

Let \(z_{i-1} = 2c_{i}\) for \(i = 1, 2, \dots, n\).

Earlier we had \(S''_i(x_{i-1}) = 2{\color{green}{c_i}}\).

So \(z_{i-1} = S''(x_{i-1})\) for \(i = 1, 2, \dots, n\). For completeness, let \(z_n = S''(x_n)\).

  • \(b_i = s_i - c_ih_i - d_ih_i^2\) for \(i = 1, 2, \dots, n\), where \(s_i = \frac{y_{i} - y_{i-1}}{h_i}\)
  • \(b_i + 2c_ih_i + 3d_ih_i^2 = b_{i+1}\) for \(i = 1, 2, \dots, n-1\)
  • \(d_i = \frac{2c_{i+1} - 2c_i}{6h_i}\) for \(i = 1, 2, \dots, n-1\)

\(d_i = \frac{z_i - z_{i-1}}{6h_i}\)

\(b_i = s_i - \frac{z_{i-1}}{2}h_i - \frac{z_i - z_{i-1}}{6h_i}h_i^2 = s_i - \left(\frac{z_{i-1}}{2} + \frac{z_i - z_{i-1}}{6}\right)h_i = s_i - \frac{z_i + 2z_{i-1}}{6}h_i\)

\(s_i - \frac{z_i + 2z_{i-1}}{6}h_i + z_{i-1}h_i + \frac{z_i - z_{i-1}}{2}h_i = s_{i+1} - \frac{z_{i+1} + 2z_i}{6}h_{i+1}\)

\(\frac{z_{i+1} + 2z_i}{6}h_{i+1} + \frac{2z_i + z_{i-1}}{6}h_i = s_{i+1} - s_i\)

Single Set of Equations

\(\left(z_{i+1} + 2z_i\right)h_{i+1} + \left(2z_i + z_{i-1}\right)h_i = 6\left(s_{i+1} - s_i\right)\) for \(n = 1, 2, \dots, n-1\).

Simplifying the left side:

\[h_iz_{i-1} + 2\left(h_i + h_{i+1}\right)z_i + h_{i+1}z_{i+1} = 6\left(s_{i+1} - s_i\right)\]

for \(n = 1, 2, \dots, n-1\).

Define \(\delta_i = h_i + h_{i+1}\) and \(B_i = 6\left(s_{i+1} - s_i\right)\) for \(i = 1, 2, \dots, n-1\).

The equations become

\[h_iz_{i-1} + 2\delta_iz_i + h_{i+1}z_{i+1} = B_i\]

Types of Splines

Two more equations:

  • Natural cubic splines: \(z_n = z_0 = 0\)

  • Periodic splines: assuming \(y_0 = y_n\), we add assumptions \(S'(x_0) = S'(x_n)\) and \(z_0 = z_n\).

  • Clamped splines: specify values for \(S'(x_0)\) and \(S'(x_n)\).

Natural Splines

The system for natural splines looks like this:

\[\begin{align*} z_0 &= 0\\ h_1z_0 + 2\delta_1z_1 + h_2z_2 &= B_1\\ h_2z_1 + 2\delta_2z_2 + h_3z_3 &= B_2\\ h_3z_2 + 2\delta_3z_3 + h_4z_4 &= B_3\\ \qquad\vdots\\ h_{n-2}z_{n-3} + 2\delta_{n-2}z_{n-2} + h_{n-1}z_{n-1} &= B_{n-2}\\ h_{n-1}z_{n-2} + 2\delta_{n-1}z_{n-1} + h_{n}z_n &= B_{n-1}\\ z_n &= 0 \end{align*}\]

Matrix form

\[ \begin{bmatrix} 1 & 0 & 0 & & & & & \\ h_1 & 2\delta_1 & h_2 & & & & & \\ &h_2 & 2\delta_2 & h_3 & & & & \\ & & \ddots & \ddots & \ddots & & & \\ & & & h_{i} & 2\delta_i & h_{i+1} & & \\ & & & & \ddots & \ddots & \ddots & \\ & & & & & h_{n-1} & 2\delta_{n-1} & h_{n} \\ & & & & & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} z_0\\z_1\\z_2\\\vdots\\z_i\\\vdots\\z_{n-1}\\z_n \end{bmatrix} = \begin{bmatrix} 0\\B_1\\B_2\\\vdots\\B_i\\\vdots\\B_{n-1}\\0 \end{bmatrix} \]

\[ \begin{bmatrix} 1 & h'_1 & 0 & & & & & \\ h_1 & 2\delta_1 & h_2 & & & & & \\ &h_2 & 2\delta_2 & h_3 & & & & \\ & & \ddots & \ddots & \ddots & & & \\ & & & h_{i} & 2\delta_i & h_{i+1} & & \\ & & & & \ddots & \ddots & \ddots & \\ & & & & & h_{n-1} & 2\delta_{n-1} & h_{n} \\ & & & & & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} z_0\\z_1\\z_2\\\vdots\\z_i\\\vdots\\z_{n-1}\\z_n \end{bmatrix} = \begin{bmatrix} \beta_0\\B_1\\B_2\\\vdots\\B_i\\\vdots\\B_{n-1}\\0 \end{bmatrix} \]

\[ \begin{bmatrix} 1 & h'_1 & 0 & & & & & \\ \phantom{h_1} & 1 & h'_2 & & & & & \\ &h_2 & 2\delta_2 & h_3 & & & & \\ & & \ddots & \ddots & \ddots & & & \\ & & & h_{i} & 2\delta_i & h_{i+1} & & \\ & & & & \ddots & \ddots & \ddots & \\ & & & & & h_{n-1} & 2\delta_{n-1} & h_{n} \\ & & & & & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} z_0\\z_1\\z_2\\\vdots\\z_i\\\vdots\\z_{n-1}\\z_n \end{bmatrix} = \begin{bmatrix} \beta_0\\\beta_1\\B_2\\\vdots\\B_i\\\vdots\\B_{n-1}\\0 \end{bmatrix} \]

\[ \begin{bmatrix} 1 & h'_1 & 0 & & & & & \\ \phantom{h_1} & 1 & h'_2 & & & & & \\ & & 1 & h'_3 & & & & \\ & & \ddots & \ddots & \ddots & & & \\ & & & h_{i} & 2\delta_i & h_{i+1} & & \\ & & & & \ddots & \ddots & \ddots & \\ & & & & & h_{n-1} & 2\delta_{n-1} & h_{n} \\ & & & & & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} z_0\\z_1\\z_2\\\vdots\\z_i\\\vdots\\z_{n-1}\\z_n \end{bmatrix} = \begin{bmatrix} \beta_0\\\beta_1\\\beta_2\\\vdots\\B_i\\\vdots\\B_{n-1}\\0 \end{bmatrix} \]

\[ \begin{bmatrix} 1 & h'_1 & & & & & & \\ & 1 & h'_2 & & & & & \\ & & 1 & h'_3 & & & & \\ & & & \ddots & \ddots & & & \\ & & & & 1 & h'_{i+1} & & \\ & & & & & \ddots & \ddots & \\ & & & & & & 1 & h'_{n} \\ & & & & & & & 1 \end{bmatrix} \begin{bmatrix} z_0\\z_1\\z_2\\\vdots\\z_i\\\vdots\\z_{n-1}\\z_n \end{bmatrix} = \begin{bmatrix} \beta_0\\\beta_1\\\beta_2\\\vdots\\\beta_i\\\vdots\\\beta_{n-1}\\\beta_n \end{bmatrix} \]

\(h'_1 = 0\), \(\beta_0 = 0\)

\(\displaystyle h'_k = \frac{h_k}{2\delta_{k-1} - h'_{k-1}h_{k-1}}\) for \(k = 2, 3, \dots, n\).

\(\displaystyle \beta_k = \frac{B_k - h_k\beta_{k-1}}{2\delta_k - h'_{k}h_{k}}\) for \(k = 1, 3, \dots, n-1\).

\(\displaystyle \beta_n = 0\)

Solve for \(z_k\)s

\[ \begin{bmatrix} 1 & h'_1 & & & & & & \\ & 1 & h'_2 & & & & & \\ & & 1 & h'_3 & & & & \\ & & & \ddots & \ddots & & & \\ & & & & 1 & h'_{i+1} & & \\ & & & & & \ddots & \ddots & \\ & & & & & & 1 & h'_{n} \\ & & & & & & & 1 \end{bmatrix} \begin{bmatrix} z_0\\z_1\\z_2\\\vdots\\z_i\\\vdots\\z_{n-1}\\z_n \end{bmatrix} = \begin{bmatrix} \beta_0\\\beta_1\\\beta_2\\\vdots\\\beta_i\\\vdots\\\beta_{n-1}\\\beta_n \end{bmatrix} \]

  • \(z_n = \beta_n = 0\)

  • \(\displaystyle z_{n-1} = \beta_{n-1} - h'_nz_n\)

  • \(\displaystyle z_{k} = \beta_{k} - h'_{k+1}z_{k+1}\)

  • \(\displaystyle z_{0} = {\color{red}{\beta_{0}}} - {\color{red}{h'_{1}}}z_{1} = {\color{red}{0}}\)

The Steps

  • Given: \(x_i\), \(y_i\), \(i = 0, 1, \dots, n\).
  • \(a_i = y_{i-1}\), \(y = 1, 2, \dots, n\).
  • \(h_i = x_i - x_{i-1}\), \(i = 1, 2, \dots, n\).
  • \(s_i = \frac{y_i - y_{i-1}}{h_i}\), \(i = 1, 2, \dots, n\).
  • \(\delta_i = h_i + h_{i+1}\), \(i = 1, 2, \dots, n-1\).
  • \(B_i = 6\left(s_{i+1} - s_i\right)\), \(i = 1, 2, \dots, n-1\).
  • \(h'_1 = 0\), \(h'_i = \frac{h_i}{2\delta_{i-1} - h'_{i-1}h_{i-1}}\), \(i = 2, 3, \dots, n\)
  • \(\beta_n = 0\), \(\beta_k = \frac{B_k - h_k\beta_{k-1}}{2\delta_k - h'_{k}h_{k}}\), \(k = n-1, n-2, \dots, 0\).
  • \(z_n = 0\), \(z_k = \beta_k - h'_{k+1}z_{k+1}\), \(k = n-1, n-2, \dots, 0\).
  • \(d_i = \frac{z_i - z_{i-1}}{6h_i}\), \(i = 1, 2,\dots, n\)
  • \(c_i = \frac{z_{i-1}}{2}\), \(i = 1, 2,\dots, n\)
  • \(b_i = s_i - \frac{z_i + 2z_{i-1}}{6}h_i\), \(i = 1, 2,\dots, n\)