Интерполяция кубическими сплайнами #

Одним из методов интерполяции сплайнами является интерполяция кубическими сплайнами.

Пусть $S(x)$ — кубический сплайн, построенный на $n$ точках $\left\{\left(x_i, y_i\right)\right\}_{i=1}^n$.
$S(x)$ состоит из $n-1$ сегмента $S_k(x)$, каждый из которых представляет из себя многочлен третьей степени.

$k$-тый сегмент кубического сплайна имеет вид:
$S_k(x) = a_k + b_k(x-x_k) + c_k(x-x_k)^2 + d_k(x-x_k)^3 {}_{, \ \forall k \in \left\{1,...,n-1\right\}}$.

Кубический сплайн принимает значения $y_i$ для каждой точки $x_i$; также он имеет непрерывные первую и вторую производные.

Общие условия кубических сплайнов: #

$S_k(x_k) = y_k {}_{, \ \forall k \in \left\{1,...,n-1\right\}}$
$S_k(x_{k+1}) = y_{k+1} {}_{, \ \forall k \in \left\{1,...,n-1\right\}}$
$S_k'(x_{k+1}) = S_{k+1}'(x_{k+1}) {}_{, \ \forall k \in \left\{1,...,n-2\right\}}$
$S_k''(x_{k+1}) = S_{k+1}''(x_{k+1}) {}_{, \ \forall k \in \left\{1,...,n-2\right\}}$

Положим шаги между узлами и значениями в узлах:

$h_k = x_{k+1} - x_k {}_{, \ \forall k \in \left\{1,...,n-1\right\}}$
$\delta_k = y_{k+1} - y_k {}_{, \ \forall k \in \left\{1,...,n-1\right\}}$.

Распишем общие условия с использованием коэффициентов:

$S_k(x_k) = y_k \Longleftrightarrow a_k = y_k$
$S_k(x_{k+1}) = y_{k+1} \Longleftrightarrow b_kh_k + c_kh_k^2 + d_kh_k^3 = \delta_k$
$S_k'(x_{k+1}) = S_{k+1}'(x_{k+1}) \Longleftrightarrow b_k + 2c_kh_k + 3d_kh_k^2 = b_{k+1}$
$S_k''(x_{k+1}) = S_{k+1}''(x_{k+1}) \Longleftrightarrow 2c_k + 6d_kh_k = 2c_{k+1}$ (можно сократить на 2)

Введём фиктивную переменную $c_n = 3d_{n-1}h_{n-1} + c_{n-1}$, что можно переписать как $c_{n-1} + 3d_{n-1}h_{n-1} = c_n$, таким образом, расширяя равенство под пунктом 4. на $k = n - 1$.

В итоге получаем:

$$ \boxed{ a_k = y_k } {}_{, \ \forall k \in \left\{1,...,n-1\right\}} \tag{1} $$$$ \boxed{ b_kh_k + c_kh_k^2 + d_kh_k^3 = \delta_k } {}_{, \ \forall k \in \left\{1,...,n-1\right\}} \tag{2} $$$$ \boxed{ b_k + 2c_kh_k + 3d_kh_k^2 = b_{k+1} } {}_{, \ \forall k \in \left\{1,...,n-2\right\}} \tag{3} $$$$ \boxed{ c_k + 3d_kh_k = c_{k+1} } {}_{, \ \forall k \in \left\{1,...,n-1\right\}} \tag{4} $$

Теперь выведем равенства, которые помогут нам выразить одни коэффициенты через другие:

$$ \boxed{ d_k = \frac{c_{k+1} - c_k}{3h_k} } {}_{, \ \forall k \in \left\{1,...,n-1\right\}} \tag{5} $$$$ \boxed{ b_k = \frac{\delta_k}{h_k} - h_k\frac{2c_k + c_{k+1}}{3} } {}_{, \ \forall k \in \left\{1,...,n-1\right\}} \tag{6} $$

Из (3), (5):
$3\frac{c_{k+1} - c_k}{3h_k}h_k^2 + 2c_kh_k + b_k = b_{k+1}$
$h_k(c_k + c_{k+1}) + b_k = b_{k+1}$

Тогда с использованием (6) можно получить следующее равенство:

Промежуточные выражения

$h_k(c_k + c_{k+1}) + \frac{\delta_k}{h_k} - h_k\frac{2c_k + c_{k+1}}3 = \frac{\delta_{k+1}}{h_{k+1}} - h_{k+1}\frac{2c_{k+1} + c_{k+2}}3$

$h_k(c_k + c_{k+1}) - h_k\frac{2c_k + c_{k+1}}3 = \frac{\delta_{k+1}}{h_{k+1}} - \frac{\delta_k}{h_k} - h_{k+1}\frac{2c_{k+1} + c_{k+2}}3$

$3h_kc_k + 3h_kc_{k+1} - 2h_kc_k - h_kc_{k+1} + 2h_{k+1}c_{k+1} + h_{k+1}c_{k+2} = 3\left(\frac{\delta_{k+1}}{h_{k+1}} - \frac{\delta_k}{h_k}\right)$

$h_kc_k + 2h_kc_{k+1} + 2h_{k+1}c_{k+1} + h_{k+1}c_{k+2} = 3\left(\frac{\delta_{k+1}}{h_{k+1}} - \frac{\delta_k}{h_k}\right)$

$$ \boxed{ c_k(h_k) + c_{k+1}(2h_k + 2h_{k+1}) + c_{k+2}(h_{k+1}) = 3\left(\frac{\delta_{k+1}}{h_{k+1}} - \frac{\delta_k}{h_k}\right) } {}_{, \ \forall k \in \left\{1,...,n-2\right\}} \tag{7} $$

Таким образом, мы получили $(n)$ переменных $c_k$ (среди них одна фиктивная — $c_n$) и $(n-2)$ уравнений для их нахождения (среди них одно получено благодаря фиктивной переменной).

На данный момент мы можем создать неполную тридиагональную матрицу, представляющую систему линейных алгебраических уравнений (СЛАУ) для нахождения значений переменных $c_k$:

$$ \begin{bmatrix} \boxed{?} & \boxed{?} & 0 & 0 & \dots & 0 & 0 \\ h_1 & 2h_1 + 2h_2 & h_2 & 0 & \dots & 0 & 0 \\ 0 & h_2 & 2h_2 + 2h_3 & h_3 & \dots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \dots & h_{n-2} & 2h_{n-2} + 2h_{n-1} & h_{n-1} \\ 0 & 0 & 0 & \dots & 0 & \boxed{?} & \boxed{?} \end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \\ c_3 \\ \vdots \\ c_{n-1} \\ c_n \end{bmatrix} = \begin{bmatrix} \boxed{?} \\ R_1 \\ R_2 \\ \vdots \\ R_{n-2} \\ \boxed{?} \end{bmatrix} \tag{8} $$

, где $R_k = 3\left(\frac{\delta_{k+1}}{h_{k+1}} - \frac{\delta_k}{h_k}\right)$.

Как видно, для формирования полной трёхдиагональной матрицы, нам необходимо ещё два уравнения.
То есть для однозначного определения коэффициентов нам необходимо два дополнительных условия.

Дополнительные условия #

I. Натуральный сплайн #

Натуральный кубический сплайн — это сплайн, у которого вторая производная равна нулю в крайних точках.
Является частным случаем “Fixed-second” сплайна.

Условия: $S_1''(x_1) = 0, \ S_{n-1}''(x_n) = 0$.

Уравнения:

$c_1 = 0$
$c_n = 0$

Матрица:

$$ \begin{bmatrix} 1 & 0 & 0 & \dots & 0 & 0 & 0 \\ h_1 & 2h_1 + 2h_2 & h_2 & \dots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \dots & h_{n-2} & 2h_{n-2} + 2h_{n-1} & h_{n-1} \\ 0 & 0 & 0 & \dots & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \\ \vdots \\ c_{n-1} \\ c_n \end{bmatrix} = \begin{bmatrix} 0 \\ R_1 \\ \vdots \\ R_{n-2} \\ 0 \end{bmatrix} $$

Подробнее про натуральный сплайн.

II. Not-a-knot сплайн #

“Not-a-knot” кубический сплайн - это сплайн, у которого третья производная непрерывна во второй и предпоследней точках. То есть третья производная равна на первом и второй сегментах, а также на предпоследнем и последнем.

Условия: $S_1'''(x_2) = S_2'''(x_2), \ S_{n-2}'''(x_{n-1}) = S_{n-1}'''(x_{n-1})$.

Уравнения:

$(h_1-h_2)c_1 + (2h_1+h_2)c_2 = 3\frac{h_1}{h_1+h_2}\left(\frac{\delta_2}{h_2} - \frac{\delta_1}{h_1}\right)$
$(h_{n-2}+2h_{n-1})c_{n-1} + (-h_{n-2}+h_{n-1})c_n = 3\frac{h_{n-1}}{h_{n-1}+h_{n-2}}\left(\frac{\delta_{n-1}}{h_{n-1}} - \frac{\delta_{n-2}}{h_{n-2}}\right)$

Матрица:

$$ \begin{bmatrix} h_1 - h_2 & 2h_1 + h_2 & 0 & \dots & 0 & 0 & 0 \\ h_1 & 2h_1 + 2h_2 & h_2 & \dots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \dots & h_{n-2} & 2h_{n-2} + 2h_{n-1} & h_{n-1} \\ 0 & 0 & 0 & \dots & 0 & h_{n-2} + 2h_{n-1} & -h_{n-2} + h_{n-1} \end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \\ \vdots \\ c_{n-1} \\ c_n \end{bmatrix} = \begin{bmatrix} R_0 \\ R_1 \\ \vdots \\ R_{n-2} \\ R_{n-1} \end{bmatrix} $$

, где $R_0 = 3\frac{h_1}{h_1+h_2}\left(\frac{\delta_2}{h_2} - \frac{\delta_1}{h_1}\right)$, $R_k = 3\left(\frac{\delta_{k+1}}{h_{k+1}} - \frac{\delta_k}{h_k}\right)$, $R_{n-1} = 3\frac{h_{n-1}}{h_{n-1}+h_{n-2}}\left(\frac{\delta_{n-1}}{h_{n-1}} - \frac{\delta_{n-2}}{h_{n-2}}\right)$.

Подробнее про “Not-a-knot” сплайн.

III. Периодический сплайн #

Периодический кубический сплайн - это сплайн, первая и вторая производные которого равна на концах.

Условия: $S_1'(x_1) = S_{n-1}'(x_n), \ S_1''(x_1) = S_{n-1}''(x_n)$.

Матрица (в данном случае не тридиагональная):

$$ \begin{bmatrix} 2h_1 + 2h_2 & h_2 & 0 & \dots & 0 & 0 & h_1 \\ h_2 & 2h_2 + 2h_3 & h_3 & \dots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \dots & h_{n-2} & 2h_{n-2} + 2h_{n-1} & h_{n-1} \\ h_1 & 0 & 0 & \dots & 0 & h_{n-1} & 2h_{n-1} + 2h_1 \end{bmatrix} \begin{bmatrix} c_2 \\ c_3 \\ \vdots \\ c_{n-1} \\ c_n \end{bmatrix} = \begin{bmatrix} R_1 \\ R_2 \\ \vdots \\ R_{n-2} \\ R_{n-1} \end{bmatrix} $$

, где $c_n = c_1$, $R_k = 3\left(\frac{\delta_{k+1}}{h_{k+1}} - \frac{\delta_k}{h_k}\right)$, $R_{n-1} = 3\left(\frac{\delta_1}{h_1} - \frac{\delta_{n-1}}{h_{n-1}}\right)$.

Подробнее про периодический сплайн.

IV. Parabolic-ends #

“Parabolic-ends” (параболические концы) кубический сплайн - это сплайн, крайние сегменты которого представляют из себя кривые второй степени или параболы.
Является частным случаем “Fixed-third” сплайна, когда третьи производные на концах устанавливаются равными нулю.

Условия: $S_1'''(x_1) = 0, \ S_{n-1}'''(x_n) = 0$.

Уравнения:

$c_1 - c_2 = 0$
$-c_{n-1} + c_n = 0$

Матрица:

$$ \begin{bmatrix} 1 & -1 & 0 & \dots & 0 & 0 & 0 \\ h_1 & 2h_1 + 2h_2 & h_2 & \dots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \dots & h_{n-2} & 2h_{n-2} + 2h_{n-1} & h_{n-1} \\ 0 & 0 & 0 & \dots & 0 & -1 & 1 \end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \\ \vdots \\ c_{n-1} \\ c_n \end{bmatrix} = \begin{bmatrix} 0 \\ R_1 \\ \vdots \\ R_{n-2} \\ 0 \end{bmatrix} $$

Особый случай: $n = 2$:

$c_1 = 0$
$c_2 = 0$

Подробнее про “Parabolic-ends” сплайн.

V. Зажатый сплайн #

Зажатый кубический сплайн - это сплайн с фиксированными первыми производными на концах.

Условия: $S_1'(x_1) = f'(x_1), \ S_{n-1}'(x_n) = f'(x_n)$.

Уравнения:

$2h_1c_1 + h_1c_2 = 3\left(\frac{\delta_1}{h_1} - f'(x_1)\right)$
$h_{n-1}c_{n-1} + 2h_{n-1}c_n = 3\left(f'(x_n) - \frac{\delta_{n-1}}{h_{n-1}}\right)$

Матрица:

$$ \begin{bmatrix} 2h_1 & h_1 & 0 & \dots & 0 & 0 & 0 \\ h_1 & 2h_1 + 2h_2 & h_2 & \dots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \dots & h_{n-2} & 2h_{n-2} + 2h_{n-1} & h_{n-1} \\ 0 & 0 & 0 & \dots & 0 & h_{n-1} & 2h_{n-1} \end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \\ \vdots \\ c_{n-1} \\ c_n \end{bmatrix} = \begin{bmatrix} R_0 \\ R_1 \\ \vdots \\ R_{n-2} \\ R_{n-1} \end{bmatrix} $$

, где $R_0 = 3\left(\frac{\delta_1}{h_1} - f'(x_1)\right)$, $R_k = 3\left(\frac{\delta_{k+1}}{h_{k+1}} - \frac{\delta_k}{h_k}\right)$, $R_{n-1} = 3\left(f'(x_n) - \frac{\delta_{n-1}}{h_{n-1}}\right)$.

Подробнее про зажатый сплайн.

VI. Fixed-second сплайн #

“Fixed-second” кубический сплайн - это сплайно с фиксированными вторыми производными на концах.
Частным случаем является натуральный сплайн, у которого вторые производные на концах устанавливаются равными нулю.

Условия: $S_1''(x_1) = f''(x_1), \ S_{n-1}''(x_n) = f''(x_n)$.

Уравнения:

$c_1 = \frac{f''(x_1)}{2}$
$c_n = \frac{f''(x_n)}{2}$

Матрица:

$$ \begin{bmatrix} 1 & 0 & 0 & \dots & 0 & 0 & 0 \\ h_1 & 2h_1 + 2h_2 & h_2 & \dots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \dots & h_{n-2} & 2h_{n-2} + 2h_{n-1} & h_{n-1} \\ 0 & 0 & 0 & \dots & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \\ \vdots \\ c_{n-1} \\ c_n \end{bmatrix} = \begin{bmatrix} f''(x_1)/2 \\ R_1 \\ \vdots \\ R_{n-2} \\ f''(x_n)/2 \end{bmatrix} $$

Подробнее про “Fixed-second” сплайн.

VII. Fixed-third сплайн #

“Fixed-third” кубический сплайн - это сплайно с фиксированными третьими производными на концах.
Частным случаем является “Parabolic-ends” сплйан, у которого третьи производные на концах устанавливаются равными нулю.

Условия: $S_1'''(x_1) = f'''(x_1), \ S_{n-1}'''(x_n) = f'''(x_n)$.

Уравнения:

$c_1 - c_2 = -\frac{h_1f'''(x_1)}{2}$
$-c_{n-1} + c_n = \frac{h_{n-1}f'''(x_n)}{2}$

Матрица:

$$ \begin{bmatrix} 1 & -1 & 0 & \dots & 0 & 0 & 0 \\ h_1 & 2h_1 + 2h_2 & h_2 & \dots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \dots & h_{n-2} & 2h_{n-2} + 2h_{n-1} & h_{n-1} \\ 0 & 0 & 0 & \dots & 0 & -1 & 1 \end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \\ \vdots \\ c_{n-1} \\ c_n \end{bmatrix} = \begin{bmatrix} -h_1f'''(x_1)/2 \\ R_1 \\ \vdots \\ R_{n-2} \\ h_{n-1}f'''(x_n)/2 \end{bmatrix} $$

Особый случай: $n = 2$:

$c_1 = -\frac{h_1\left(f'''(x_1) + f'''(x_2)\right)}{8}$
$c_2 = -c_1$

Подробнее про “Fixed-third” сплайн.

Произвольные условия #

Мы можем установить два произвольных условия для двух заданных точек.

В этом случае мы можем построить “подсплайн”, начальной и конечной точками которого будут эти две точки.
Примечание: могут быть тонкости для Fixed-third условий, для которых указываются сегменты, и для Not-a-knot условий, для которых указываются вторая и предпоследняя точки.

В этих точках можно задать необходимые уравнения, соответствующие условиям для первой или последней точки (см. вышеописанные варианты).

Далее мы можем итеративно построить предыдущие и последующие сегменты, расширяя подсплайн до концов сплайна.

Ниже представлены формулы итеративного расчёта коэффициентов $k$-го сегмента при условии известных коэффициентов $(k+1)$-го или $(k-1)$-го.

Формулы для итеративного построения #

Даны: $a_{k+1}, b_{k+1}, c_{k+1}, d_{k+1}$. Вывести $a_k, b_k, c_k, d_k$. $_{k \in \left\{1,...,n-2\right\}}$.

$d_k = \frac{c_{k+1} - \frac{b_{k+1}}{h_k} + \frac{\delta_k}{h_k^2}}{h_k}$.
$c_k = c_{k+1} - 3d_kh_k$.
$b_k = \frac{\delta_k}{h_k} - d_kh_k^2 - c_kh_k$.
$a_k = y_k$.

Даны: $a_{k-1}, b_{k-1}, c_{k-1}$. Вывести $a_k, b_k, c_k$. $_{k \in \left\{2,...,n-1\right\}}$.

$a_k = y_k$.
$b_k = 3d_{k-1}h_{k-1}^2 + 2c_{k-1}h_{k-1} + b_{k-1}$.
$c_k = 3d_{k-1}h_{k-1} + c_{k-1}$.
$d_k = \frac{\frac{\delta_k}{h_k^2} - c_k - \frac{b_k}{h_k}}{h_k}$.

Примечания #

1. Условия “Зажатый” (Clamped) и “Fixed-second” в одной точке

Общая схема вычислений не покрывает случай, когда в одной точке заданы два условия: “Зажатый” и “Fixed-second”.

Выведем коэффициенты отдельно для этого случая.

Пусть в точке $x_k$ заданы первая производная $f'$ и вторая производная $f''$.

Коэффициенты предыдущего ($k-1$-го) сегмента:

$a_{k-1} = y_{k-1}$
$b_{k-1} = \frac{3\delta_{k-1}}{h_{k-1}} + \frac{f''h_{k-1}}{2} - 2f'$
$c_{k-1} = \frac{3f'}{h_{k-1}} - \frac{3\delta_{k-1}}{h_{k-1}^2} - f''$
$d_{k-1} = \frac{f''}{2h_{k-1}} - \frac{f'}{h_{k-1}^2} + \frac{\delta_{k-1}}{h_{k-1}^3}$

Коэффициенты текущего ($k$-го) сегмента:

$a_k = y_k$
$b_k = f'$
$c_k = \frac{f''}{2}$
$d_k = \frac{\delta_k}{h_k^3} - \frac{f'}{h_k^2} - \frac{f''}{2h_k}$

2. Условия “Fixed-second” и “Not-a-knot” в одной точке

Коэффициенты для $k-1$-го и $k$-го сегментов:

$a_{k-1} = y_{k-1}$
$a_k = y_k$
$c_k = \frac{f''(x_k)}{2}$
$c_{k-1} = \frac{f'' - 6dh_{k-1}}{2}$
$b_{k-1} = \frac{\delta_{k-1}}{h_{k-1}} - \frac{f''h_{k-1}}{2} + 2dh_{k-1}^2$
$b_k = \frac{\delta_{k-1}}{h_{k-1}} + \frac{f''h_{k-1}}{2} - dh_{k-1}^2$
$d_{k-1} = d_k = d = \frac{\frac{\delta_k}{h_k} - \frac{\delta_{k-1}}{h_{k-1}} - \frac{f''(x_k)}{2}h_{k-1} - \frac{f''(x_k)}{2}h_k}{h_k^2 - h_{k-1}^2}$

Как видно, для вычисления коэффициентов нам необходимо выполнить деление на выражение $h_k^2 - h_{k-1}^2$, которое равно нулю, если шаги $h_{k-1}$ и $h_k$ равны.

Таким образом, в частности, при интерполяции с использованием равномерно распределённых точек, установить оба условия “Fixed-second” и “Not-a-knot” в одной точке невозможно.