Cubic interpolation
For cubic interpolation, the polynomial order is raised to \(n=3\). Now, for a
given interal \(\bar x \in [0, 1]^N\) we're looking for the values
\(a_{i_1,\dots,i_N}\) in the following function
\[
s(\bar x_1,\dots,\bar x_N) = \sum_{i_1,\dots,i_N=0}^3 \bar a_{i_1 \dots i_N}\prod_{k=0}^N\bar x_{k}^{i_k}
\]
Expressing \(s(\bar x_1,\dots,\bar x_N)\) in terms of \(f(\bar x_1,\dots,\bar x_N)\)
leads to:
\[
\begin{equation}
s(\bar x_1,\dots,\bar x_N) = \sum_{i_1,\dots,i_N=0}^1 \sum_{j_1,\dots,j_N=0}^1 f^{(j_1,\dots,j_N)}(i_1,\dots,i_N) \prod_{k=1}^N \alpha^{(3,j_k)}_{i_k} (\bar x_k), \label{eq:s_normalized}
\end{equation}
\]
where
\[
\begin{align*}
\alpha_0^{3,0}(\bar x_k) &= 1 - 3\bar x_k^2 + 2 \bar x_k ^3\\
\alpha_0^{3,1}(\bar x_k) &= 3\bar x_k^2 - 2\bar x_k^3\\
\alpha_1^{3,0}(\bar x_k) &= \bar x_k - 2 \bar x_k^2 + \bar x_k^3\\
\alpha_1^{3,1}(\bar x_k) &=-\bar x_k^2 + \bar x_k^3
\end{align*}
\]
The reader is reminded that, as stated here, \(f^{(l)}(\bar x) \equiv
\frac{\textrm{d}^lf(\bar x)}{\textrm{d}\bar x^l}\). For cubic interpolation only
zeroth and first order differentation occurs, which occasionaly might be denoted
using \(f(\bar x)\) and \(f'(\bar x)\) respectively.
Analogue to linear interpolation, \(\eqref{eq:s_normalized}\) is also rewritten to its non-normalized version:
\[
s(x_1,\dots,x_N) = \sum_{i_1,\dots,i_N=0}^3 a_{i_1 \dots i_N}\prod_{k=0}^N x_{k}^{i_k}
\]
in which we're looking for expressions for the coefficients \(a_{i_1 \dots i_N}\).
By substututing \(\bar x_k = (x_k - {}^0x_k)/h_k\) (where \(h_k={}^1x_k-{}^0x_k\))
into \(\eqref{eq:s_normalized}\) we obtain
\[
\begin{equation}
\begin{split}
s(x_1,\dots,x_N) = \sum_{i_1,\dots,i_N=0}^1 \sum_{j_1,\dots,j_N=0}^1 &
\left( \prod_{k=1}^N h_k^{j_k} \right) f^{(j_1,\dots,j_N)}({}^{i_1}x_1,\dots,{}^{i_N}x_N) \dots \\
& \dots \prod_{k=1}^N \alpha^{(3,j_k)}_{i_k} (x_k)
\end{split} \label{eq:s_non_normalized}
\end{equation}
\]
Note the additional term \(\prod_{k=1}^N h_k^{j_k}\), which arises due to
substition of \(\bar x_k = (x_k - {}^0x_k)/h_k\) into \(f^{(l)}(\bar x_k)\).
In \(\eqref{eq:s_non_normalized}\) \(\alpha_{i_k}^{(3,j_k)}\) can be expressed as:
\[
\begin{align*}
\alpha_0^{3,0}(x_k) &= \frac{1}{h_k^3} \sum^3_{i=0} \delta^{(0,i)}_k x_k^i\\
\alpha_0^{3,1}(x_k) &= \frac{1}{h_k^3} \sum^3_{i=0} \delta^{(1,i)}_k x_k^i\\
\alpha_1^{3,0}(x_k) &= \frac{1}{h_k^3} \sum^3_{i=0} \delta^{(2,i)}_k x_k^i\\
\alpha_1^{3,1}(x_k) &= \frac{1}{h_k^3} \sum^3_{i=0} \delta^{(3,i)}_k x_k^i
\end{align*}
\]
in which the expressions \(\delta^{(j,i)}_k\) read
|
\(i=0\) |
\(i=1\) |
\(i=2\) |
\(i=3\) |
| \(j=0\) |
\({}^1x_k^2({}^1x_k-3{}^0x_k)\) |
\(+6{}^0x_k{}^1x_k\) |
\(-3({}^0x_k+{}^1x_k)\) |
\(+2\) |
| \(j=1\) |
\(-{}^0x_k{}^1x_k^2\) |
\({}^1x_k(2{}^0x_k+{}^1x_k)\) |
\(-({}^0x_k+2{}^1x_k)\) |
\(+1\) |
| \(j=2\) |
\({}^0x_k^2(3{}^1x_k-{}^0x_k)\) |
\(-6{}^0x_k{}^1x_k\) |
\(+3({}^0x_k+{}^1x_k)\) |
\(-2\) |
| \(j=3\) |
\(-{}^1x_k{}^0x_k^2\) |
\({}^0x_k({}^0x_k+2{}^1x_k)\) |
\(-(2{}^0x_k+{}^1x_k)\) |
\(+1\) |
Note that expression for \(\alpha_i^{3,j}\) can be written even shorter by introducing a binary number \(ij\), i.e. \(00, 01, 10, 11 = 0, 1, 2, 3\):
\[
\begin{equation}
\alpha_i^{3,j}(x_k) = \frac{1}{h_k^3} \sum_{l=0}^3 \delta_k^{(ij,l)}x_k^i \label{eq:unified_alpha}
\end{equation}
\]
The following sections show how to employ \(\eqref{eq:s_non_normalized}\) to
obtain expressions for the coefficients \(a_{i_1 \dots i_N}\) for 1, 2 and N
dimensions.
1 dimension
For 1 dimension the interval definition simplifies to \(x \in [x_0, x_1]\). With
this:
\[
\begin{equation}
s(x) = \sum_{k=0}^3 a_k x^k = a_0 + a_1 x + a_2 x^2 + a_3 x^3
\end{equation}
\]
in which the coefficients can be written as
\[
\begin{align*}
a_0 &= (f(x_ 0)x_1^2(x_1 - 3 x_0) + f(x_1)x_0^2(3x_1 - x_0) - h x_0 x_1(x_1f'(x_0) + x_0f'(x_1)))/h^3\\
a_1 &= (+6 x_0 x_1 (f_0-f_1) + h ( x_1 (2 x_0 + x_1)f'(x_0) + x_0 (x_0 + 2 x_1)f'(x_1)))/h^3\\
a_2 &= (-3 (x_0 + x_1)(f_0-f_1) - h( (x_0 + 2 x_1)f'(x_0) + (2 x_0 + x_1)f'(x_1)))/h^3\\
a_3 &= (+2 (f_0-f_1) + h (f'(x_0) + f'(x_1)))/h^3
\end{align*}
\]
where \(h = x_1 - x_0\).
2 dimensions
For 2 dimensions the interval definition now reads \(x \in [{}^0x, {}^1x]^2\).
With this:
\[
\begin{align*}
s(x_1,x_2) = & \sum_{k,l=0}^3 a_{kl} x_1^k x_2^l
\end{align*}
\]
Now, we're looking for the coefficients \(a_{kl}\) for \(k,l=0,1,2,3\).
Analogue to linear interpolationt this is accomplished by rewriting
\(\eqref{eq:s_non_normalized}\) for \(N=2\). To simplify the expressions, we shall
rewrite \(f_{ij}^{(k,l)} = f^{(k,l)}({}^{i}x_1{}^{j}x_2)\). This leads
to
\[
\begin{align*}
s(x_1,x_2) = & \sum_{i_1,i_2=0}^1 \sum_{j_1,j_2=0}^1
h_1^{j_1}h_2^{j_2} f_{i_1i_2}^{(j_1,j_2)} \alpha^{(3,j_1)}_{i_1} (x_1) \alpha^{(3,j_2)}_{i_2} (x_2)\\
=& \frac{1}{h_1^3 h_2^3} \sum_{i_1,i_2=0}^1 \sum_{j_1,j_2=0}^1
h_1^{j_1}h_2^{j_2} f_{i_1i_2}^{(j_1,j_2)} \sum_{k,l=0}^3\delta_1^{(i_1j_1,k)} \delta_2^{(i_2j_2,l)} \\
=& \frac{1}{h_1^3 h_2^3} \sum^3_{k,l=0} \Bigl( \\
& f_{00}^{(0,0)} \delta^{(0,k)}_1 \delta^{(0,l)}_2 +
f_{01}^{(0,0)} \delta^{(0,k)}_1 \delta^{(1,l)}_2 +
f_{10}^{(0,0)} \delta^{(1,k)}_1 \delta^{(0,l)}_2 +
f_{11}^{(0,0)} \delta^{(1,k)}_1 \delta^{(1,l)}_2 ~ + \\
& h_2\left(f_{00}^{(0,1)} \delta^{(0,k)}_1 \delta^{(2,l)}_2 +
f_{01}^{(0,1)} \delta^{(0,k)}_1 \delta^{(3,l)}_2 +
f_{10}^{(0,1)} \delta^{(1,k)}_1 \delta^{(2,l)}_2 +
f_{11}^{(0,1)} \delta^{(1,k)}_1 \delta^{(3,l)}_2\right) ~ + \\
& h_1\left(f_{00}^{(1,0)} \delta^{(2,k)}_1 \delta^{(0,l)}_2 +
f_{01}^{(1,0)} \delta^{(2,k)}_1 \delta^{(1,l)}_2 +
f_{10}^{(1,0)} \delta^{(3,k)}_1 \delta^{(0,l)}_2 +
f_{11}^{(1,0)} \delta^{(3,k)}_1 \delta^{(1,l)}_2 \right)~ + \\
& h_1h_2\left(f_{00}^{(1,1)} \delta^{(2,k)}_1 \delta^{(2,l)}_2 +
f_{01}^{(1,1)} \delta^{(2,k)}_1 \delta^{(3,l)}_2 +
f_{10}^{(1,1)} \delta^{(3,k)}_1 \delta^{(2,l)}_2 +
f_{11}^{(1,1)} \delta^{(3,k)}_1 \delta^{(3,l)}_2 \right)\\
& \Bigr)x_1^k x_2^l
\end{align*}
\]
From this it becomes clear that
\[
\begin{align*}
a_{kl} = & \frac{1}{h_1^3 h_2^3}
\Biggl( f_{00}^{(0,0)} \delta^{(0,k)}_1 \delta^{(0,l)}_2 +
f_{01}^{(0,0)} \delta^{(0,k)}_1 \delta^{(1,l)}_2 +
f_{10}^{(0,0)} \delta^{(1,k)}_1 \delta^{(0,l)}_2 +
f_{11}^{(0,0)} \delta^{(1,k)}_1 \delta^{(1,l)}_2 ~ + \\
& h_2\left(f_{00}^{(0,1)} \delta^{(0,k)}_1 \delta^{(2,l)}_2 +
f_{01}^{(0,1)} \delta^{(0,k)}_1 \delta^{(3,l)}_2 +
f_{10}^{(0,1)} \delta^{(1,k)}_1 \delta^{(2,l)}_2 +
f_{11}^{(0,1)} \delta^{(1,k)}_1 \delta^{(3,l)}_2\right) ~ + \\
& h_1\left(f_{00}^{(1,0)} \delta^{(2,k)}_1 \delta^{(0,l)}_2 +
f_{01}^{(1,0)} \delta^{(2,k)}_1 \delta^{(1,l)}_2 +
f_{10}^{(1,0)} \delta^{(3,k)}_1 \delta^{(0,l)}_2 +
f_{11}^{(1,0)} \delta^{(3,k)}_1 \delta^{(1,l)}_2 \right)~ + \\
& h_1h_2\left(f_{00}^{(1,1)} \delta^{(2,k)}_1 \delta^{(2,l)}_2 +
f_{01}^{(1,1)} \delta^{(2,k)}_1 \delta^{(3,l)}_2 +
f_{10}^{(1,1)} \delta^{(3,k)}_1 \delta^{(2,l)}_2 +
f_{11}^{(1,1)} \delta^{(3,k)}_1 \delta^{(3,l)}_2 \right)\Biggr)\\
\end{align*}
\]
N dimensions
For N dimensions, replacing the \(\alpha\)-terms in
\(\eqref{eq:s_non_normalized}\) with the simplified \(\alpha\)-terms from
\(\eqref{eq:unified_alpha}\):
\[
\begin{equation}
\begin{split}
s(x_1,\dots,x_N) = \sum_{i_1,\dots,i_N=0}^1 \sum_{j_1,\dots,j_N=0}^1 &
\left( \prod_{k=1}^N h_k^{j_k} \right) f^{(j_1,\dots,j_N)}({}^{i_1}x_1,\dots,{}^{i_N}x_N) \dots \\
& \dots \prod_{k=1}^N \frac{1}{h_k^3} \sum_{l=0}^3 \delta_k^{(i_kj_k,l)}x_k^i
\end{split}
\end{equation}
\]
This can be reordered to
\[
s(x_1,\dots,x_N) = \sum_{\mathbf{m}\in\{0,1,2,3\}^N} A_{\textbf{m}}\prod_{k=1}^N
x_k^{m_k},
\]
where
\[
A_{\textbf{m}} =\sum_{i_1,\dots,i_N=0}^1 \sum_{j_1,\dots,j_N=0}^1
f^{(j_1,\dots,j_N)}\Bigl({}^{i_1}x_1,\dots,{}^{i_N}x_N\Bigr)
\prod_{k=1}^N \frac{h_k^{\,j_k}}{h_k^3}\,\delta_k^{(i_kj_k,m_k)}
\]
In here, \(\textbf{m}\) is defined as
\[
\mathbf{m}=(m_1,\dots,m_N) \quad \text{with} \quad m_k\in\{0,1,2,3\}.
\]
\(A_\textbf{m}\) can be rewritten to the original coefficients: \(A_\textbf{m} =
A_{(m_1,\dots,m_N)} = a_{m_1 \dots m_N}\).
Role of Derivatives in Cubic Interpolation
The formulas in the preceding sections show that cubic polynomial coefficients depend on both function values \(f({}^{i}x_1,\dots,{}^{i}x_N)\) and first derivatives \(f'({}^{i}x_1,\dots,{}^{i}x_N)\) at grid points. For example, in the 1D case the coefficients \(a_0, a_1, a_2, a_3\) explicitly depend on derivatives \(f'(x_0)\) and \(f'(x_1)\).
While function values are given from the input data, first derivatives at grid points must be computed. The Derivatives section provides several methods for estimating these derivatives, each with different properties and characteristics.