Интерполяция сплайнами

Пусть задана таблица значений функции y_i в узлах х₀ < х₁ < ... < х_п.Обозначим h_i = x_i – x_i-1, i = 1, 2, ... , п.

Сплайн– гладкая кривая, проходящая через задан­ные точки (х_i, y_i), i = 0, 1, ... , п. Интерполяция сплай­нами заключается в том, что на каждом отрезке [х_i-1, x_i]используется многочлен определенной степени. Наиболее часто применяется многочлен третьей степени, реже – второй или четвертой. При этом для определения коэффициентов многочленов используются условия непре­рывности производных в узлах интерполяции.

Интерполяция кубическими сплайнамипредставляет собой локальную интерполяцию, когда на каждом отрез­ке [х_i-1, x_i], i = 1, 2, ... , п применяется кубическая кри­вая, удовлетворяющая некоторым условиям гладкости, а именно, непрерывности самой функции и ее первой и вто­рой производных в узловых точках. Использование куби­ческой функции вызвано следующими соображениями. Если предположить, что интерполяционная кривая соот­ветствует упругой линейке, закрепленной в точках (х_i, y_i),то из курса сопротивления материалов известно, что эта кривая определяется как решение дифференциального уравнения f^(IV)(x) = 0 на отрезке [х_i-1, x_i](для простоты из­ложения мы не рассматриваем вопросы, связанные с физи­ческими размерностями). Общим решением такого уравне­ния является многочлен 3-й степени с произвольными коэффициентами, который удобно записать в виде
S_i(x) = а_i + b_i(х - x_i-1) + с_i(x - x_i-1)²+ d_i(x - x_i-1)³,
х_i-1 £ х £ х_i, i = 1, 2, ... , п.(4.32)

Коэффициенты функций S_i(x)определяются из усло­вий непрерывности функции и ее первой и второй произ­водных во внутренних узлах x_i, i = 1, 2,..., п - 1.

Из формул (4.32) при х = х_i-1 получим

S_i(x_i-1) = y_i-1 = a_i, i = 1, 2,..., п,(4.33)

а при х = х_i

S_i(x_i) = а_i + b_ih_i + с_ih_i²+ d_ih_i³,(4.34)

i = 1, 2,..., n.

Условия непрерывности интерполяционной функции записываются в виде S_i(x_i) = S_i-1(x_i), i = 1, 2, ... , n - 1 и из условий (4.33) и (4.34) следует, что они выполнимы.

Найдем производные функции S_i(x):

S'_i(x) = b_i + 2с_i(х - x_i-1) + 3di(х – x_i-1)²,

S"_i(x) = 2c_i + 6d_i(x - x_i-1).

При x = x_i-1, имеем S'_i(x_i-1) = b_i, S" (x_i-1) = 2с_i, а при х = х_i получим

S'_i(x_i) = b_i + 2с_ih_i + 3dih_i², S" (x_i) = 2с_i + 6d_ih_i.

Условия непрерывности производных приводят к уравнениям

S'_i(x_i) = S'_i+1(x_i) Þ b_i + 2с_ih_i + 3dih_i² = b_i+1,

i = l, 2,... , п - 1. (4.35)

S"_i (x_i) = S"_i+1(x_i) Þ 2с_i + 6d_ih_i = 2c_i+1,

i = l, 2,..., n - 1. (4.36)

Всего имеем 4n – 2 уравнений для определения 4n не­известных. Чтобы получить еще два уравнения, исполь­зуют дополнительные краевые условия, например, требо­вание нулевой кривизны интерполяционной кривой в концевых точках, т. е. равенства нулю второй производ­ной на концах отрезка [а, b] а = х₀, b = х_n:

S"₁(x₀) = 2c₁ = 0 Þ с₁ = 0,

S"_n(x_n) = 2с_n + 6d_nh_n = 0 Þ с_n + 3d_nh_n = 0. (4.37)

Систему уравнений (4.33)–(4.37) можно упростить и получить рекуррентные формулы для вычисления коэф­фициентов сплайна.

Из условия (4.33) имеем явные формулы для вычисле­ния коэффициентов a_i:

a_i = y_i-1, i= 1,..., n. (4.38)

Выразим d_i через c_i с помощью (4.36), (4.37):

; i = 1, 2,...,n; .

Положим с_n+1= 0, тогда для d_i получим одну формулу:

, i = 1, 2,...,n. (4.39)

Подставим выражения для а_i и d_i в равенство (4.34):

, i = 1, 2,..., n.

и выразим b_i, через с_i:

, i = 1, 2,..., n. (4.40)

Исключим из уравнений (4.35) коэффициенты b_i и d_i с помощью (4.39) и (4.40):

,

i = 1, 2,..., n -1.

Отсюда получим систему уравнений для определения с_i:

(4.41)

Система уравнений (4.41) может быть переписана в виде

(4.42)

Здесь введено обозначение

, i =1, 2,..., n - 1.

Решим систему уравнений (4.42) методом прогонки. Из первого уравнения выразим с₂через с₃:

c₂ = a₂c₃ + b₂, , . (4.43)

Подставим (4.43) во второе уравнение (4.42):

h₂(a₂c₃ + b₂) + 2(h₂ + h₃)c₃ + h₃c₄ = g₂,

и выразим с₃ через с₄:

с₃ = a₃с₄ + b₃, (4.44)

Предполагая, что с_i-1= a_i-1c_i + b_i-1 из i-го уравне­ния (4.42) получим

c_i = a_iс_i+1 + b_i

, i = 3,..., n – 1, a_n = 0, (4.45)

, i = 3,..., n.

Сформулируем алгоритм интерполяции с помощью кубического сплайна.

Исходные данные: значения функции у₀, у₁,..., у_п в узлах х₀, х₁,..., х_п (х₀ < х₁< ... < х_n).

0. Вычислим значения h_i и g_i:

h_i = x_i – x_i-1 i = 1, 2,..., n,

, i =1, 2,..., n - 1. (4.46)

1. Прямой ход прогонки для вычисления коэффици­ентов c_i. Вычислим коэффициенты прогонки a_i и b_i:

, ,

, i = 3,..., n – 1, a_n = 0, (4.47)

, i = 3,..., n.

2. Обратный ход прогонки для вычисления коэффици­ентов с_i. Вычислим коэффициенты c_i:

c_n+1 = 0,

c_i = a_iс_i+1 + b_i, i = n, n -1,..., 2, (4.48)

c₁ = 0.

3. Вычисление коэффициентов а_i, b_i, d_i:

a_i = y_i-1,

, (4.49)

i = 1, 2,..., n.

4. Вычисление значения функции с помощью сплай­на. Для этого найти такое значение i,что данное значе­ние переменной х принадлежит отрезку [x_i-1, x_i] и вы­числить

S_i(x) = а_i + b_i(х - x_i-1) + с_i(x - x_i-1)²+ d_i(x - x_i-1)³. (4.50)