Пусть каждому значению аргумента xi, i= 0 ,...,n соответствуют значения функции f (xi) =yi и требуется найти функциональную зависимость в виде сплайна (2.4.1), удовлетворяющего перечисленным ниже требованиям:
1) функция S 3(x) непрерывна на отрезке [ a, b ] вместе со своими производными до второго порядка включительно;
2) S 3(xi)= yi, i= 0,1 ,...,n;
3) функция S 3(x) удовлетворяет одному из вариантов краевых условий а-в (условия а,б,в могут задаваться смешанно, т.е. на одном конце первая производная, на другом - вторая).
Сформулированная задача имеет единственное решение.
Для понижения порядка системы уравнений, которую приходится решать для определения коэффициентов сплайна, кубические многочлены записываются в специальном виде, что позволяет использовать большую часть условий при построении.
Вторая производная S" 3(x) непрерывна и на каждом отрезке [ xi -1, xi ], (i= 1,..., n), выражается линейной функцией, поэтому представим ее в виде интерполяционного многочлена Лагранжа 1-й степени
, (2.4.4)
где hi = xi - xi -1, mi = S" 3(xi).
Проинтегрируем обе части равенства (2.4.4)
, (2.4.5)
где 
. В силу непрерывности 
.
Отсюда
, i= 1, ...,n (2.4.6)
Проинтегрируем (2.4.5)
,
(2.4.7)
где 
. В силу непрерывности 
,.
Тем самым, получаем систему уравнений
, i= 1, ...,n. (2.4.8)
Последовательно вычитая последующее уравнение из предыдущего, получим
,
или с учетом (2.4.6)
, i= 1, ...,n -1. (2.4.9)
Если известны m 0= f" (a) и mn = f" (b), то mi, i= 1, ...,n -1 определяются из решения системы линейных алгебраических уравнений (2.4.9).
Пусть известны m 0= f" (a) и pn = f¢ (b). Тогда систему (2.4.9) необходимо дополнить. Согласно (2.4.8) n -е уравнение
. (2.4.10)
Если известны mn = f" (b) и p 0= f¢ (a), то подставляя полученное из (2.4.6) значение p 1 в (2.4.8) при i =1 найдем
|
|
.
Тем самым получается уравнение с номером i =0, дополняющее систему (2.4.9)
. (2.4.11)
Если известны p 0= f¢ (a) и pn = f¢ (b), то к системе (2.4.9) добавляются оба уравнения (2.4.10) и (2.4.11). В случае задания краевых условий вида в) к системе (2.4.9) вместо (2.4.10) и (2.4.11) добавляются уравнения, полученные дифференцированием (2.4.4)
, (2.4.12)
Система линейных алгебраических уравнений (2.4.9)-(2.4.12) имеет трехдиагональную матрицу с диагональным преобладанием. Такие матрицы являются неособенными. Поэтому неизвестные m 0, m 2 ,..., mn находятся однозначно.
В результате решения системы уравнений определяются параметры mi, i =0,…, n. Тогда выразив pi из (2.4.8) и подставив в (2.4.7), получим
(2.4.13)
По этой формуле вычисляются значения функции S 3(x).
Решение системы линейных алгебраических уравнений может быть найдено итерационными и прямыми методами решения, в том числе методом прогонки. Для решения задачи этим методом требуется выполнить примерно 3 n сложений, 3 n умножений, 2 n делений.
Метод прогонки
Пусть имеется система уравнений, записанная в матричном виде:
. (2.4.14)
При использовании краевых условий а (на обеих границах заданы значения первой производной) согласно (2.4.9)-(2.4.11)
,
, 
,
.
При использовании краевых условий типа б и в изменяются только значения ai, bi, ci с индексами i =0 и i = n. Если на обеих границах заданы значения второй производной (краевые условия б)
,
.
При использовании краевых условий в (на обеих границах заданы значения третьей производной) согласно (2.4.12)
,
.
Решение системы ищется в виде
mi = l i mi +1 + m i, i =0,..., n -1, (2.4.15)
где l i, m i - прогоночные коэффициенты. Используя выражение для mi- 1 из (2.4.15), подставим это неизвестное в i ‑е уравнение системы
|
|
cimi- 1 + aimi + bimi+ 1 = di.
Получаем
(ai +ci l i- 1) mi + bi mi+ 1 = di -ci m i- 1, 
.
Сравнивая это соотношение с (2.4.15), выводим рекуррентные формулы для прогоночных коэффициентов l i, m i (прямая прогонка):
.
(2.4.16)
Очевидно, что mn= m n (при bn= 0). Все остальные неизвестные находим по формулам (2.4.15), используя выражения для прогоночных коэффициентов (2.4.16).
Численный эксперимент
В первых четырех колонках табл. 2.5.1 представлены результаты интерполяции функции sin x на отрезке [0,p] при равномерном разбиении с удвоением числа отрезков n. Краевые условия – равенство нулю второй производной (б). В колонке D max представлена максимальная погрешность | S 3(x)- f (x)|, вычисленная в точках, находящихся между узлами сетки; K D - отношение погрешности предыдущей строки к данной (коэффициент уменьшения погрешности при удвоении n). Видно, что K D сохраняет значение, соответствующее четвертому порядку точности (K D»24) до значений n =1000-3000, выше которых в общей погрешности результата преобладает погрешность округления.
Табл. 2.5.1
| n | Dmax | Dоц | K D | n | Dmax | Dоц | K D | |
| 4.47×10-04 | - | - | 1.07×10-01 | - | - | |||
| 2.57×10-05 | 2.80×10-05 | 17.42 | 5.33×10-02 | 5.00×10-02 | 2.01 | |||
| 1.59×10-06 | 1.60×10-06 | 16.15 | 2.67×10-02 | 2.49×10-02 | 2.00 | |||
| 9.92×10-08 | 9.94×10-08 | 16.04 | 1.33×10-02 | 1.24×10-02 | 2.00 | |||
| 6.19×10-09 | 6.20×10-09 | 16.01 | 6.67×10-03 | 6.22×10-03 | 2.00 | |||
| 3.87×10-10 | 3.87×10-10 | 16.00 | 3.33×10-03 | 3.11×10-03 | 2.00 | |||
| 2.42×10-11 | 2.42×10-11 | 16.00 | 1.67×10-03 | 1.56×10-03 | 2.00 | |||
| 1.51×10-12 | 1.51×10-12 | 16.00 | 8.33×10-04 | 7.78×10-04 | 2.00 | |||
| 9.46×10-14 | 9.45×10-14 | 15.99 | 4.17×10-04 | 3.89×10-04 | 2.00 | |||
| 6.00×10-15 | 5.91×10-15 | 15.77 | 2.08×10-04 | 1.95×10-04 | 2.00 | |||
| 4.97×10-16 | 3.75×10-16 | 12.07 | 1.04×10-04 | 9.73×10-05 | 2.00 | |||
| 1.77×10-16 | 2.78×10-17 | 2.82 | 5.21×10-05 | 4.86×10-05 | 2.00 |
|
|
Отметим, что для практического применения оценки погрешности (2.4.3) необходимо знать верхнюю оценку k -й производной функции f (x). Это не всегда возможно, так как сплайн часто используется для представления искомой функции, которая является решением дифференциального или интегрального уравнения. Для оценки погрешности можно применить правило, использующее закономерность зависимости погрешности от h или n. Однако необходимо заметить, что при увеличении числа узлов погрешность интерполяции в какой-нибудь конкретной точке x может изменяться нерегулярно, поскольку положение этой точки относительно соседних узлов (отношение (x - xj -1)/(xj - xj -1)) может различаться для разных n. Как показывает численный эксперимент, максимальная погрешность на отрезке [ a, b ] уменьшается в K D»2 k раз при удвоении n. Используя свойство сохранения значения K D для максимума погрешности Dmax(n) можно получить оценку в виде
. (2.5.1)
Для этого необходимо иметь способ оценки величины Dmax(n), если даже точное значение интерполируемой функции f (x) неизвестно. Можно использовать следующий способ, заключающийся в сравнении значений S 3(x), вычисленных при разных числах отрезков, на которые разбивается отрезок [ a, b ] - n и 2 n. При удвоении n при равномерном или неравномерном разбиении возникает n новых узловых точек xj -1/2, лежащих между общими узлами xj -1 и xj (рис. 2.5.1). Тогда в качестве оценки Dmax(n) выберем
. (2.5.2)
Как видно из табл. 2.5.1, величина Dоц(n) с достаточной точностью оценивает погрешность Dmax(n) для каждого n при значениях n £2520, (выше которых преобладает погрешность округления).
Рис. 2.5.1. К оценке погрешности при интерполяции сплайном
В последних четырех колонках табл. 2.5.1 представлены результаты аналогичных расчетов для интерполяции той же функции с краевыми условиями 
. Как показывают результаты, ошибка в задании краевых значений (в данном случае точные значения производной 
) приводит к погрешности интерполяции, порядок которой (в данном случае первый) совпадает с порядком погрешности, обусловленной разрывом соответствующей производной.