Сплайн-функции. Интерполиране с кубични сплайни. Теорема на Холидей.

Малко предговор, за да си имаме на предства що е това сплайн-фунцкия.
Идеята се ражда преди време, когато машините (компютрите) са имали проблеми с пресмятането на полиноми от висoка степен. Затова учените са се сетили, че могат да интерполират дадена функция с полином от сравнително ниска степен, стига интервалът, в който функцията се интерполира, да е достатъчно малък. Когато се изиска интерполацията да става върху голям интервал, той се разбива на по-малки (достатъчно малки) интервали и за всеки отделен интервал се намира полином (който е от ниска степен), интерполиращ функцията.

Дефиниция (Сплайн-функция)

Нека $x_1 < x_2 < \dots < x_n$ са различни реални числа.
Сплайн-функция от степен $r$ с възли $x_1, x_2, \dots, x_n$ наричаме функцията $s(x)$, притежаваща следните свойства

(1)
\begin{eqnarray} &(i)& s(x)|_{x \in (x_i, x_{i+1})} \in \pi_r, \ i = 0, 1, 2, \dots, n,\ x_0 = -\infty,\ x_{n+1} = \infty \\ &(ii)& s(x) \in C^{r-1}(-\infty, \infty) \end{eqnarray}

Идеята на първото свойство е получената полиномиална крива да е "съставена" от полиноми, чиято степен не надвишава $r$.
Колкото до второто условие - неговата задача е постигането на гладкост в крайните точки на малките интервали.

С $S_r(x_1, x_2, \dots, x_n)$ ще отбелязваме класа от всички сплайн-функции от степен $r$ с възли $x_1, x_2, \dots, x_n$.

Дефиниция (отсечена степенна функция)

Отсечената степенна функция е сплайн-функция с един възел.
Дефинира се по следния начин:

(2)
\begin{align} (x - \xi)^r_{+} := \begin{cases} (x - \xi)^r, & x \geq \xi \\ 0, & x < \xi \end{cases} \end{align}

Пример при $r = 3, \xi = 0$:
base_spline

Теорема

Всяка сплайн-функция $s(x) \in S_r(x_1, x_2, \dots, x_n)$ може да се представи по единствен начин във вида

(3)
\begin{align} s(x) = p_0(x) + \sum_{k = 1}^n c_k(x - x_k)^r_{+} \end{align}

където $p_0 \in \pi_r$ и числата $\{c_k\}_{k = 1}^n$ се получават от формулата

(4)
\begin{align} c_k = \frac{s^{(r)}(x_k + 0) - s^{(r)}(x_k - 0)}{r!} \end{align}

Доказателство:
Нека $s(x) \in S_r(x_1, x_2, \dots, x_n)$
Тогава $s(x) \equiv p_i(x) \in \pi_r$ при $x \in (x_i, x_{i+1})$ за $i = 0, 1, \dots, n$. (Където $p_i$ е полином)
В частност, при $x \in (x_0, x_1) = (-\infty, x_1),\ s(x) = p_0(x) \in \pi_r$, понеже $(x - x_k)^r_{+} = 0, \ k = 1, 2, \dots, n$

Имаме, че
$s(x) = p_{i-1}(x), \ x \in (x_{i-1}, x_i)$
$s(x) = p_i(x),\ x \in (x_i, x_{i+1})$

Освен това, обаче, $s(x)$ има непрекъсната $(r-1)$-ва производна в интервала $(-\infty, \infty)$.
Тогава е изпълнено, че

(5)
\begin{align} p^{(j)}_{i-1}(x_i) = p^{(j)}_i(x_i), \ j = 0, 1, \dots, r-1 \end{align}

$\Longrightarrow p_i(x) - p_{i-1}(x)$ има $r$ -кратна нула в точката $x_i$,т.е
$p_i(x) - p_{i-1}(x) = c_i(x - x_i)^r \Longrightarrow p_i(x) = p_{i-1}(x) + c_i(x - x_i)^r_{+} \quad x \in (x_i, x_{i+1})$, където $c_i$ е някакъв коефициент, разбира се.

Минавайки надясно аналогично се получава и формулата:

(6)
\begin{align} s(x) = p_0(x) + \sum_{i = 1}^n c_i(x - x_i)^r_{+} \end{align}

имайки предвид, че

(7)
\begin{align} c_i(x - x_i)^r_{+} = \begin{cases} c_i(x - x_i)^r, & x \in (x_i, x_{i+1}) \\ 0, & x \in (x_{i-1}, x_i) \end{cases} \end{align}

Сега остава да покажем, че коефициентите са еднозначно определими.
Имаме, че

(8)
\begin{equation} p_i(x) - p_{i-1}(x) = c_i(x - x_i)^r \end{equation}

Диференцираме $r$ пъти и вземаме стойността в точката $x_i$ получаваме

(9)
\begin{equation} p_i^{(r)}(x_i) - p_{i-1}^{(r)}(x_i) = c_i r! \end{equation}

и като си спомним, че всъщност $s(x) \equiv p_{i-1}(x) \And s(x) \equiv p_i(x)$ в съответните интервали, то:

(10)
\begin{align} s^{(r)}(x_i + 0) - s^{(r)}(x_i - 0) = c_i r! \Longrightarrow c_i = \frac{s^{(r)}(x_i + 0) - s^{(r)}(x_i - 0)}{r!} \end{align}

с което теоремата е доказана. $\Box$

Следствие

Като следствие от теоремата можем да твърдим, че функциите $1, x, x^2, \dots, x^r, (x-x_1)^r_{+}, (x-x_2)^r_{+}, \dots, (x-x_n)^r_{+}$ са базис на $S_r(x_1, x_2, \dots, x_n)$. Откъдето следва, че размерността му е точно броят на тези функции, а именно: $n + r + 1$.

На части кубична интерполация

Тази точка няма нищо общо със сплайни. Като за начало, това което ще получим няма да има непрекъсната втора производна, което е задължително условие за сплайн от 3та степен (ака - кубичен). От друга страна начинът по който извеждаме формулата много прилича на същинското интерполиране със кубичен сплайн в следващата точка. Идеята на това тук повече прилича на съставните квадратурни формули на правоъгълниците/трапците/Симпсон - просто разбираме интервала на части и във всяка част опитваме да приближим функцията с полином. В тази точка ще приближаваме всяка една част със полином от 3та степен, като ще гарантираме съвпадение на нулевата и първата производна на функцията, с нейната интерполираща функция (не е нито полином, нито е сплайн, а е комбинация от полиноми подобно на сплайн).

Нека $f(x)$ е дефинирана върху крайния затворен интевал $[a,b]$.
Търсим полином $s(x)$, такъв, че $s(x)|_{(x_i, x_{i+1})} = p_i \in \pi_3$ и, освен това, да е изпълнено, че:

(11)
\begin{array} {|cccccc} p_i(x_i) & = &f(x_i) & =: & f_i \\ p'_i(x_i) & = & f'(x_i) & =: & f'_i \\ p_i(x_{i+1}) & = & f(x_{i+1}) & =: & f_{i+1} \\ p'_i(x_{i+1}) & = & f'(x_{i+1}) & =: & f'_{i+1} \end{array}

Това е една типична интерполационна задача на Ермит.
Сега, за да намерим всеки полином $p_i(x)$, използваме интерполационната формула на Нютон, която използва разделени разлики с кратни възли и изглежда ето така:

(12)
\begin{equation} p_i(x) = p_i[x_i] + p_i[x_i, x_i](x - x_i) + p_i[x_i, x_i, x_{i+1}] (x - x_i)^2 + p_i[x_i, x_i, x_{i+1}, x_{i + 1}] (x - x_i)^2 (x - x_{i+1}) \end{equation}

Определяме си разделените разлики от дефиницията както следва:

(13)
\begin{eqnarray} p[x_i] & = & p(x_i) = f(x_i) = f_i \\ p[x_i, x_i] & = & \frac{p'(x_i)}{1!} = f'(x_i) = f'_i \\ p[x_i, x_i, x_{i+1}] & = & \frac{p[x_i, x_{i+1}] - p[x_i, x_i]}{x_{i+1} - x_i} = \frac{\frac{f_{i+1} - f_i}{x_{i+1} - x_i} - f'_i}{x_{i+1} - x_i}\\ p[x_i, x_i, x_{i+1}, x_{i+1}] & = & \frac{p[x_i, x_{i+1}, x_{i+1}] - p[x_i, x_i, x_{i+1}]}{ x_{i+1} - x_i} = \dots = \frac{f'_{i+1} - 2 \frac{f_{i+1} - f_i}{x_{i+1} - x_i} + f'_i}{(x_{i+1} - x_i)^2} \end{eqnarray}

Така вече нямаме неизвестни във формулата за $p_i(x)$ и сме супер щастливи.

Оценка на грешката

При $x \in [x_i, x_{i+1}]$ имаме, че грешката при Ермит се смята както при интерполацията по Лагранж
$f(x) - s(x) = f(x) - p_i(x) = \frac{f^{(4)}(\xi)}{4!}(x - x_i)^2 (x - x_{i+1})^2$

Правим си следните означения за удобство:
$\Delta_i := x_{i+1} - x_i$
$\Delta := \max \Delta_i, \ 0 \leq i \leq n -1$

Ако $|f^{(4)}(x) | \leq M$ в $[a, b]$, тогава е изпълнено, че

(14)
\begin{align} | f(x) - s(x)| \leq \frac{M}{24} \underbrace{\max_{x \in [x_i, x_{i+1}]} (x - x_i)^2 (x - x_{i+1})^2}_{\leq \Delta^4/16} \leq \frac{M}{384} \Delta^4 \end{align}

Ако интервалът $[a,b]$ е разделен на равни подинтервали, т.е $x_{i+1} - x_i = \frac{b - a}{n}, \ i = 0, 1, \dots, n-1$, тогава ще имаме

(15)
\begin{align} |f(x) - s(x)| \leq \frac{M(b-a)^4}{384n^4} \end{align}

което клони към 0 при $n \to \infty$, както $\frac{1}{n^4}$

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

Нека имаме точките $a = x_0 < x_1 < \dots <x_{n-1} < x_n = b$. и функция $f(x)$, която е достатъчно гладка в интервала $[a, b]$.
Търсим кубичен слайн $s(x)|_{(x_i, x_{i+1})} = p_i (x) \in \pi_3$, който интерполира функцията в тези точки.
Т.е изпълнени са следните условия за $i = 0, 1, \dots, n-1$:

(16)
\begin{array} {|cccccc} p_i(x_i) & = & f_i\\ p'_i(x_i) & = & s_i\\ p_i(x_{i+1}) & = & f_{i+1} \\ p'_i(x_{i+1}) & = & s_{i+1} \end{array}

Тук $s_i$ и $s_{i+1}$ са все още неизвестни, но после ще ги определим, така че да е изпълнено сплайн функцията да има непрекъснати производни до 2ра степен включително. Т.е. това което правим е да параметризираме стойностите на първата производна (чрез $s_i$) за да може И втората производна да се получи равна. Ако фиксираме някакви конкретни стойности ще стане като в предишната точка.

Тъй като горе-написаната система е интерполационна задача на Ермит (а ние знаем какво прави Ермит в такъв случай… нали?), остава да си сметнем разделените разликивъв възлите $x_i, x_i, x_{i+1}$ и $x_{i+1}$, за да получим коефициентите на $p_i(x)$1

И така, полиномът $p_i(x)$ има следния вид:

(17)
\begin{align} p_i(x) = f_i + s_i(x - x_i) + \bigg(\frac{\frac{f_{i+1} - f_i}{\Delta_i} - s_i}{\Delta_i}\bigg) (x - x_i)^2 + \frac{s_{i+1} - 2 \frac{f_{i+1} - f_i}{\Delta_i} + s_i}{\Delta_i^2} (x - x_i)^2(x - x_{i+1}) \end{align}

Тук, разбира се, $\Delta_i = x_{i+1} - x_i$.

Сега остава да удовлетворим и условията за гладкост:

(18)
\begin{align} p_{i - 1}''(x_i) = p_i''(x_i), \ i = 0, 1, \dots, n-1 \end{align}

Тъй като имаме полиномите в явен вид, можем директно да сметнем вторите им производни.

(19)
\begin{eqnarray} p_i''(x_i) &=& 2 \frac{\frac{f_{i+1} - f_i}{\Delta_i} - s_i}{\Delta_i} - 2 \frac{s_{i+1} - 2 \frac{f_{i+1} - f_i}{\Delta_i} + s_i}{\Delta_i} \\ p_{i-1}''(x_i) &=& 2 \frac{\frac{f_i - f_{i-1}}{\Delta_{i-1}} - s_{i-1}}{\Delta_{i-1}} + 4 \frac{s_i - 2\frac{f_i - f_{i-1}}{\Delta_{i-1}} + s_{i-1}}{\Delta_{i-1}} \end{eqnarray}

..който не вярва - да си ги сметне сам ;)

Сега ги приравняваме, прехвърляме неизвестните $s_0, s_1, \dots, s_n$ от лявата страна и получаваме следната система

(20)
\begin{eqnarray} & &\frac{1}{\Delta_{i-1}}s_{i-1} + \bigg(\frac{2}{\Delta_{i-1}} + \frac{2}{\Delta_i}\bigg)s_i + \frac{1}{\Delta_i}s_{i+1} = 3 \bigg( \frac{f_i - f_{i-1}}{\Delta_{i-1}^2} + \frac{f_{i+1} - f_i}{\Delta_i^2} \bigg) \\ & &i = 1, 2, \dots, n-1 \end{eqnarray}

Матрицата на тази система е тридиагонална, т.е има ненулеви елементи само по главния диагонал и диагоналите под и над него. Както след малко ще покажем, тази матрица е с доминиращ главен диагонал и винаги има решение, което се намира по метода на прогонката, което също ще видим по-късно (в следващата тема).


Видове кубични сплайни

Както видяхте по-горе кубичен сплайн се получава, като решим една тридиагонална матрица, от която намираме първите производни във възлите, от където използвайки Ермит намираме полиномите от 3та степен за всеки интервал. Параметрите които избирахме за първата производна - $s_i$ ще бъдат фиксирани от решението на матрицата с изключение на $s_0, s_n$. Сега ще покажем два начина, по които могат да се изберат и ще дадем имена на съответстващите им сплайни

Естествен кубичен сплайн

Ако подберем така $s_0, s_n$, че $s''(a) = s''(b) = 0$ (т.е сплайна в крайните точки има нулева втора производна) то така избраният кубичен сплайн се нарича естествен кубичен сплайн.

Пълна кубична сплайн интерполация

Ако функцията която интерполираме има непрекъсната първа производна: $f \in C^1[a, b]$, то може да подберем $s_0 = f'(a) \quad s_n = f'(b)$. Така избраният кубичен сплайн се нарича пълна кубична сплайн интерполация.

Теорема (на Холидей)

Нека $\bar{x} = (x_0, x_1, \dots, x_n), \ x_0 = a, \ x_n = b$ и $\bar{y} = (y_0, y_1, \dots, y_n)$ и нека с $F(\bar{x}, \bar{y})$ означим всички функции $f \in C^2[a,b]$, такива че $f(x_i) = y_i, \ i = 0, 1, \dots, n$

Тогава за всяка функция $f(x) \in F(\bar{x}, \bar{y})$ е изпълнено, че

(21)
\begin{align} \int_a^b [f''(x)]^2 \, dx \geq \int_a^b [s''(x)]^2 \, dx \end{align}

където $s(x)$ е естественият кубичен сплайн от $F(\bar{x}, \bar{y})$.

Равенството се постига т.с.т.к $f(x) \equiv s(x)$

Доказателство:
Като за начало да спомена за същестуването на следната

Лема на Питагор:

Ако $\int_{a}^{b} g(x)h(x) dx = 0$, тогава $\int_{a}^{b} [g(x) + h(x)]^2 dx \geq \int_{a}^{b} g^2(x) dx$
Доказателство:

(22)
\begin{eqnarray} \int_a^b [g(x) + h(x)]^2\ dx &=& \int_a^b [g^2(x) + 2g(x)h(x) + h^2(x)]\ dx \\ &=& \int_a^b g^2(x)\ dx + 2\underbrace{\int_a^b g(x)h(x)\ dx}_{= 0} + \underbrace{\int h^2(x)\ dx}_{\ge 0} \\ &\ge& \int_a^b g^2(x)\ dx \\ \end{eqnarray}

Ако $h(x)$ е непрекъсната, то равенство се получава при $h(x) \equiv 0$ в интервала $[a, b]$.

Сега да си продължим доказателството на теоремата. Ще представим $f''(x) = f''(x) - s''(x) + s''(x)$ и ще приложим лемата на Питагор за функциите $f''(x) - s''(x)$ и $s''(x)$:

(23)
\begin{align} \int_a^b [f''(x)]^2\ dx = \int_a^b [f''(x) - s''(x) + s''(x)]^2\ dx \ge \int_a^b [s''(x)]^2\ dx \end{align}

Разбира се, не сме доказали най-важното - а именно че $\int_a^b [f''(x) - s''(x)]s''(x)\ dx = 0$ (това беше необходимо условие за да приложим Лемата).

(24)
\begin{eqnarray} \int_a^b [f''(x) - s''(x)]s''(x)\ dx &=& \sum_{k=0}^{n-1} \int_{x_k}^{x_{k+1}} [f''(x) - s''(x)]s''(x)\ dx \\ &=& \sum_{k=0}^{n-1} \int_{x_k}^{x_{k+1}} s''(x)\ d[f'(x) - s'(x)] \\ &=& \sum_{k=0}^{n-1} s''(x)[f'(x) - s'(x)]\Big|_{x_k}^{x_{k+1}} - \sum_{k=0}^{n-1} \int_{x_k}^{x_{k+1}}[f'(x) - s'(x)]\underbrace{s'''(x)}_{c_k} dx \\ &=& \underbrace{s''(b)}_{=0}[f'(b) - s'(b)] - \underbrace{s''(a)}_{=0}[f'(a) - s'(a)] - \sum_{k=0}^{n-1} c_k \int_{x_k}^{x_{k+1}}[f'(x) - s'(x)]\ dx \\ &=& - \sum_{k=0}^{n-1} c_k [f(x) - s(x)]\Big|_{x_k}^{x_{k+1}} \\ &=& 0 \\ \end{eqnarray}

Използвахме, че $s''(a) = s''(b) = 0$ (това е от дефиницията на естествен кубичен сплайн), $s'''(x) = \mathrm{const}$ за $x \in (x_i, x_{i+1})$, защото във всеки един отделен интервал функцията е всъщност полином от 3та степен, чиито 3та производна е константа (6 * коефициента пред $x^3$). Накрая ползвахме и $f(x) = s(x)$ във всички възли.
Равенство може да има само при $f''(x) - s''(x) = 0$ т.с.т.к. $f'(x) - s'(x) = \mathrm{const}$ т.с.т.к. $f(x) - s(x) \in \pi_1$, но функциите са равни в точките $a, b$ следователно ако и двете бяха полиноми от първа степен щяха да съвпадат. Следователно равенство само при $f(x) \equiv s(x)$.$\Box$

Забележка: Неравенството от теоремата: $\int_a^b [f''(x)]^2\ dx \ge \int_a^b[s''(x)]^2\ dx$ е валидно и за класа

(25)
\begin{align} \tilde F(\bar x, \bar y) = \{ f \in C^2[a, b] : f(x_i) = y_i \quad i = \overline{0,n}, f'(a) = A, f'(b) = B \} \end{align}

Където $A, B$ са дадени константи и $\tilde s(x)$ е пълната кубична сплайн интерполация (просто там имаме, че $f'(a) - s'(a) = f'(b) - s'(b) = 0$ и 2те събираеми пак се унищожават). Обърнете внимание, че на този клас (както и на горния) отговаря точно един сплайн.

Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-NonCommercial-ShareAlike 3.0 License