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

Пусть задана таблица значений функции y_i в узлах x ₀ < x ₁ < … < x_n. Обозначим h_i = x_i – x_{i – 1}, i = 1, 2, …, n.

Сплайн (spline) — гладкая кривая, проходящая через заданные точки (x_i, y_i), i = 0, 1, …, n. Интерполяция сплайнами представляет собой локальную интерполяцию, когда на каждом отрезке [ x_{i – 1}, x_i ], i = 1, 2, …, n, применяется кубическая кривая, удовлетворяющая некоторым условиям гладкости, а именно, непрерывности самой функции и её первой и второй производных в узловых точках. Использование кубической функции вызвано следующими соображениями. Если предположить, что интерполяционная кривая соответствует упругой линейке, закрепленной в точках (x_i, y_i), то из курса сопротивления материалов известно, что эта кривая определяется как решение дифференциального уравнения f ^(IV)(x) = 0 на отрезке [ x_{i – 1}, x_i ] (для простоты изложения мы здесь не рассматриваем вопросы, связанные с физическими размерностями). Общим решением такого уравнения является многочлен 3-й степени с произвольными коэффициентами, который удобно записать в виде

S_i (x) = a_i + b_i (x – x_{i – 1}) + c_i (x – x_{i – 1})² + d_i (x – x_{i – 1})³, (4.32)

x_{i – 1} ≤ x ≤ x_i, i = 1, 2, …, n.

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

Из (4.32) при x = x_{i – 1} получим

S_i (x_{i – 1}) = y_{i – 1} = a_i, i = 1, 2, …, n, (4.33)

а при x = x_i

S_i (x_i) = a_i + b_ih_i + c_ih_i ² + d_ih_i ³ = y_i, i = 1, 2, …, n. (4.34)

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

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

При x = x_{i – 1} имеем , а при x = x_i получим

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

(4.35)

(4.36)

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

(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):

Положим c_{n + 1} = 0 и тогда для d_i получим одну формулу

(4.39)

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

и выразим b_i через c_i:

(4.40)

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

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

(4.41)

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

(4.42)

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

.

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

(4.43)

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

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

(4.44)

Предполагая, что из i -го уравнения (4.42) получим

(4.45)

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

Исходные данные:

Значения функции y ₀, y ₁, …, y_n в узлах x ₀, x ₁, …, x_n (x ₀ < x ₁ < … < x_n).

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

h_i = x_i – x_{i – 1}, i = 1, 2, …, n,

. (4.46)

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

(4.47)

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

(4.48)

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

(4.49)

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

S_i (x) = a_i + b_i (x – x_{i – 1}) + c_i (x – x_{i – 1})² + d_i (x – x_{i – 1})³. (4.50)

Пример 4.6. Построить кубический сплайн для функции f (x) = sin(π x) на отрезке [0; 2], используя разбиение отрезка на n = 10 частей. Найти значение в точке x = 0,48.

Решение. Применим для решения задачи в программе Excel формулы (4.42) — (4.46). В диапазоне A 1: M 1 запишем обозначения столбцов как в таблице 4.6.

0. Построим таблицу значений функции.

В ячейках A 2: A 12 запишем значения индекса i = 0, 1, …, 10.

Шаг изменения x равен 2/10 = 0,2. В ячейку B 2 запишем 0, в B 3 — число 0,2; выделим B 2: B 3 и маркером заполнения протянем вниз до B 12. Присвоим диапазону B 2: B 12 имя x командой «Вставка — Имя — Присвоить — x».

В ячейку C 3 запишем формулу =B3–B2; выделим C 3 и маркером заполнения протянем вниз до C 12.

В ячейку D 2 запишем формулу =SIN(3,1415926*x); выделим D 2 и маркером заполнения протянем вниз до D 12.

В ячейку E 3 запишем формулу =3*(D4-D3)/C4-3*(D3-D2)/C3; выделим E 3 и маркером заполнения протянем вниз до E 11.

1. Вычислим коэффициенты прогонки.

В ячейку F 4 запишем формулу =-C4/(2*(C3+C4)); В ячейку F 5 запишем формулу =(-C5-F4*C4)/(2*(C4+C5)); выделим F 5 и маркером заполнения протянем вниз до F 11; В ячейку F 12 запишем 0.

В ячейку G 4 запишем формулу =E3/(2*(C4+C5)); В ячейку G 5 запишем формулу =(E4-G4*C4)/(2*(C5+C6)). Выделим G 5 и маркером заполнения протянем вниз до F 12.

2. Вычислим коэффициенты c_i.

В ячейку H 13 запишем число 0. В ячейку H 12 запишем формулу =F12*H13+G12. Выделим H 12 и маркером заполнения протянем вверх до H 4. В ячейку H 3 запишем число 0.

3. Вычислим коэффициенты сплайна.

В ячейку I 3 запишем формулу =D2. Выделим I 3 и маркером заполнения протянем вниз до I 12.

В ячейку J 3 запишем формулу =(D3-D2)/C3-(2*H3+H4)*C3/3. Выделим J 3 и маркером заполнения протянем вниз до J 12.

В ячейку K 3 запишем формулу =(H4-H3)/(3*C3). Выделим K 3 и маркером заполнения протянем вниз до K 12.

4. Вычисление значения сплайна.

Точка x = 0,48 попадает в отрезок [0,4; 0,6]. Следовательно, надо использовать S ₃(x), т.е. строку i = 3. Поэтому в ячейку L 5 запишем формулу =0,48-0,4. В ячейку M 5 запишем формулу =I5+J5*L5+H5*L5^2+K5*L5^3.

Получим значение 0,997. Точное значение sin(0,48π) = 0,998.

Табл. 4.6

	A	B	C	D	E	F	G	H	I	J	K	L	M
	i	x	h	y	g	α	β	c	a	b	d	x-xi-1	S(x)
		0,000		0,000
		0,200	0,200	0,588	-3,37			0,00	0,00	3,14	-4,99
		0,400	0,200	0,951	-5,45	-0,25	-4,21	-2,99	0,59	2,54	-3,11
		0,600	0,200	0,951	-5,45	-0,19	-5,76	-4,86	0,95	0,97	0,13	0,08	0,997
		0,800	0,200	0,588	-3,37	-0,20	-5,37	-4,78	0,95	-0,99	3,16
		1,000	0,200	0,000	0,00	-0,20	-2,87	-2,89	0,59	-2,56	4,99
		1,200	0,200	-0,588	3,37	-0,20	0,72	0,11	0,00	-3,16	4,88
		1,400	0,200	-0,951	5,45	-0,20	4,03	3,04	-0,59	-2,55	3,20
		1,600	0,200	-0,951	5,45	-0,20	5,80	4,96	-0,95	-0,94	-1,25
		1,800	0,200	-0,588	3,37	-0,20	5,36	4,21	-0,95	0,87	2,54
		2,000	0,200	0,000		0,00	5,74	5,74	-0,59	2,17	-9,57
								0,00

Дата добавления: 2015-07-20