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

Лекция 3. Интерполяция сплайнами. Пусть дана табличная функция: i 1 ... n xi x0 x1 . . . xn yi y0 y1 . . . yn или yi = f (xi ), i = 0, 1, . . . , n. (1) Интерполяционные формулы Ньютона при использовании большого числа узлов интерполяции на всем отрезке [x0 , xn ] часто приводят к плохому приближению из-за накопления погрешностей в процессе вычислений. Для снижения погрешностей отрезок [x0 , xn ] разбивается на частичные отрезки и на каждом из них функцию заменяют приближенно полиномом невысокой степени. Сплайном называется кусочно-полиномиальная непрерывная функция, определенная на отрезке [x0 , xn ] и имеющая на этом отрезке некоторое количество непрерывных производных. Сплайны стали широко использоваться в вычислительной математике сравнительно недавно. В машиностроительном черчении они применяются уже давно, так как сплайны – это лекала или гибкие линейки, деформация которых позволяет провести кривую через заданные точки (xi , yi ), i = 0, 1, . . . , n. Чаще всего применяется интерполирование функции кубическим сплайном. Кубический сплайн, соответствующий узлам интерполяции (1), – функция S(x), удовлетворяющая следующим условиям: 1) на каждом отрезке [xi−1 , xi ], i = 1, 2, . . . , n, функция S(x) является кубическим многочленом; 2) функция S(x), а также ее первая и вторая производные непрерывны на отрезке [x0 , xn ]; 3) S(xi ) = yi , i = 0, 1, 2, . . . , n (условия интерполяции). Рассмотрим способ построения кубического сплайна. 1 На каждом из отрезков [xi−1 , xi ], i = 1, 2, . . . , n, будем искать сплайнфункцию S(x) = Si (x) в виде полинома третьей степени: Si (x) = ai + bi (x − xi ) + ci (x − xi )2 + di (x − xi )3 , x ∈ [xi−1 , xi ], где ai , bi , ci , di , i = 1, 2, . . . , n, – искомые коэффициенты (рис. 1). y 6 s Si−1 (x) s s S1 (x) s Sn (x) s s O x0 s Si (x) s x1 s xi−1 xi xn - x Рис. 1: Тогда Si0 (x) = bi + 2ci (x − xi ) + 3di (x − xi )2 , Si00 (x) = 2ci + 6di (x − xi ), откуда Si (xi ) = ai , Si0 (xi ) = bi , Si00 (xi ) = 2ci , i = 1, 2, . . . , n (2) Условия непрерывности функции и ее первой и второй производных записываются в виде Si (xi−1 ) = Si−1 (xi−1 ), Si0 (xi−1 ) = Si−1 (xi−1 ), 00 Si00 (xi−1 ) = Si−1 (xi−1 ), i = 2, . . . , n, (3) а условия интерполяции в виде Si (xi ) = yi , i = 1, 2, . . . , n, S1 (x0 ) = y0 . (4) Таким образом, имеем 3(n − 1) + n + 1 = 4n − 2 уравнения для определения 4n коэффициентов сплайна ai , bi , ci , di , i = 1, 2, . . . , n. Для однозначного задания сплайна перечисленных условий недостаточно. Необходимо наложить два дополнительных условия. Например, граничные условия: S 00 (x0 ) = S100 (x0 ) = 0, S 00 (xn ) = Sn00 (xn ) = 0 2 (5) (такой сплайн называется естественным). Пусть hi = xi − xi−1 , i = 1, 2, . . . , n, a0 = y0 . Тогда, используя (2), запишем систему уравнений (3) – (5) для определения коэффициентов сплайна в следующем виде: ai = y i , i = 1, 2, . . . , n; (6) ai − bi hi + ci h2i − di h3i = ai−1 , bi − 2ci hi + 3di h2i = bi−1 , 2ci − 6di hi = 2ci−1 , 2c1 − 6d1 h1 = 0, i = 1, . . . , n; i = 2, . . . , n; i = 2, . . . , n; (7) (8) (9) 2cn = 0. (10) Данную систему линейных алгебраических уравнений можно решить любым известным методом. Однако в целях экономии памяти ЭВМ и машинного времени эту систему можно привести к более удобному виду. Пусть c0 = 0. Тогда из (9), (10) выражаем di : di = ci − ci−1 , 3hi i = 1, . . . , n. (11) Подставляя (6) и (11) в (7), находим коэффициенты bi : bi = yi − yi−1 2ci + ci−1 + hi , hi 3 i = 1, . . . , n. (12) С помощью (11) и (12) исключаем из (8) коэффициенты bi и di : ci−2 hi−1 + 2ci−1 (hi + hi−1 ) + ci hi = y − y yi−1 − yi−2  i i−1 =3 − , hi hi−1 c0 = cn = 0. i = 2, . . . , n, (13) Таким образом, получен следующий алгоритм построения кубического сплайна. Алгоритм построения кубического сплайна: 1) Находим ai = yi , i = 1, 2, . . . , n. 3 2) Вычисляем hi = xi − xi−1 , i = 1, 2, . . . , n, и находим ci из системы линейных уравнений: ci−1 hi + 2ci (hi+1 + hi ) + ci+1 hi+1 = y − y yi − yi−1  i+1 i , − =3 hi+1 hi c0 = cn = 0. i = 1, 2, . . . , n − 1, ci − ci−1 , i = 1, 2, . . . , n. 3hi yi − yi−1 2ci + ci−1 + hi , i = 1, 2, . . . , n. 4) Находим bi = hi 3 5) Записываем ответ для сплайна в виде: 3) Находим di = S1 (x) = a1 + b1 (x − x1 ) + c1 (x − x1 )2 + d1 (x − x1 )3 , x ∈ [x0 , x1 ], ... Sn (x) = an + bn (x − xn ) + cn (x − xn )2 + dn (x − xn )3 , x ∈ [xn−1 , xn ]. Пример. Построить кубический сплайн для функции y(x), заданной таблично: i 0 1 2 3 xi 1 3 5 7 yi 4 -2 6 -3  Будем искать сплайн-функцию S(x) в виде: S1 (x) = a1 + b1 (x − 3) + c1 (x − 3)2 + d1 (x − 3)3 , x ∈ [1, 3], S2 (x) = a2 + b2 (x − 5) + c2 (x − 5)2 + d2 (x − 5)3 , x ∈ [3, 5], S3 (x) = a3 + b3 (x − 7) + c3 (x − 7)2 + d3 (x − 7)3 , x ∈ [5, 7]. 1) Находим ai : a1 = y1 = −2, a2 = y2 = 6, a3 = y3 = −3. 2) Имеем hi = xi − xi−1 = 2, i = 1, 2, 3, и система уравнений для коэффициентов ci имеет вид:   y − y y − y 1 2 1   − , c0 h1 + 2c1 (h2 + h1 ) + c2 h2 = 3   h h  2 1 y − y y 2 − y1  3 2 − , c1 h2 + 2c2 (h3 + h2 ) + c3 h3 = 3  h h  3 2    c0 = 0, c3 = 0 4 или     2c0 + 8c1 + 2c2 = 21, ( 2c1 + 8c2 + 2c3 = − 51 2 ,    c = 0, 8c1 + 2c2 = 21, 2c1 + 8c2 = − 51 2 , c3 = 0 откуда c2 = −4,1. c1 = 3,65, 3) Находим di : d1 = c1 − c0 ≈ 0,61, 3h1 d2 = c2 − c1 ≈ −1,30, 3h2 d3 = c3 − c2 ≈ 0,68. 3h3 4) Находим bi : b1 = y1 − y0 2c1 + c0 + h1 ≈ 1,87, h1 3 b2 = y2 − y1 2c2 + c1 + h2 ≈ 0,97, h2 3 y3 − y2 2c3 + c2 + h3 ≈ −7,23. h3 3 5) Искомый кубический сплайн имеет вид:   −2 + 1,87(x − 3) + 3,65(x − 3)2 + 0,61(x − 3)3 ,    S(x) = 6 + 0,97(x − 5) − 4,1(x − 5)2 − 1,30(x − 5)3 ,     −3 − 7,23(x − 7) + 0,68(x − 7)3 , b3 = 5 x ∈ [1, 3], x ∈ [3, 5], x ∈ [5, 7].