4.2. Интерполяционный многочлен Ньютона

Имеем случай неравностоящих узлов, n = 3;

N₃(x) = f(x₀) + (x–x₀)f(x₀,x₁) + (x–x₀)(x–x₁)f(x₀,x₁,x₂) + (x–x₀)(x–x₁)(x–x₂)f(x₀,x₁,x₂,x₃).

По схеме таблицы 2 находим раздельные разности

f(x₀,x₁) = ;

f(x₁,x₂) = ;

f(x₂,x₃) = ;

f(x₀,x₁,x₂) =

f(x₁,x₂,x₃) =

f(x₀,x₁,x₂,x₃) = .

Результаты расчетов поместим в таблицу:

n	xn	fn	f(xn, xn+1)	f(xn, xn+1, xn+2)	f(xn, xn+1, xn+2, xn+3)
0	0	–0,5
1	0,1	0	5	–40/3	125/3
2	0,3	0,2	1	15/2
3	0,5	1	4

Используя первые в столбцах разделенные разности, получим

N₃(x) = –0,5 + (x – 0)5 + (x – 0)(x – 0,1)(– ) + (x – 0)(x – 0,1)(x – 0,3) =

= x³ – 30x² + x – 0,5 . (30)

Аналогично расчету по Лагранжу.

Напомним, что расчеты интерполяционного многочлена Ньютона выполняются по формуле

,

где – текущая точка, в которой надо вычислить значение многочлена;

– разделенные разности порядка k, которые вычисляются по следующим рекуррентным формулам:

Схема алгоритма расчета многочлена Ньютона, реализованная в виде функции PN с параметрами, значения которых аналогичны рассмотренной ранее функции PL, представлена на рис. 5.2.

Результатом функции PN является значение N.

Рис. 5.2. Схема расчета многочлена Ньютона

5.5. Сплайны

Пусть интервал [a,b] разбит узлами x_i, как и выше, на n отрезков, 0  i  n Сплайном S_n(x) называется функция, определенная на [a, b], принадлежащая C^k[a, b] и такая, что на каждом отрезке [x_i, x_i+1], 0  i  n–1 – это полином n-й степени.

В частности, это могут быть, построенные специальным образом, многочлены 3-й степени (кубический сплайн), которые являются математической моделью гибкого тонкого стержня, закрепленного в двух точках на концах с заданными углами наклона  и .

В данной физической модели стержень принимает форму, минимизирующую его потенциальную энергию. Пусть форма стержня определяется какой-то функцией y = S(x). Из курса сопротивления материалов известно, что уравнение свободного равновесия имеет вид S^(IV)(x) = 0. А этому состоянию соответствует многочлен третьей степени между двумя соседними узлами интерполяции. Его выбирают в виде

S(x) = a_i + b_i(x – x_i–1) + c_i(x – x_i–1)² + d_i(x – x_i–1)³ ; x_i–1  х  x_i . (31)

Стоит проблема нахождения a_i, b_i, c_i, d_i. Для определения их на всех n элементарных участках интервала [a,b] необходимо составить 4n уравнений. Часть этих уравнений в составе 2n получают из условия прохождения S(x) через заданные точки, т.е.

S(x_i–1) = y_i–1; S(x_i) = y_i .

Эти условия можно записать, используя (31) в виде:

(32)

(33)

Уравнения в количестве (2n–2) получают из условия непрерывности первых и вторых производных в узлах интерполяции. Условие гладкости.

Вычислим производные многочлена (31)

(x) = b_i + 2c_i(x – x_i–1) + 3d_i(x – x_i–1)²,

(x) = 2c_i + 6d_i(x – x_i–1); при x_i–1  х  x_i . (34)

Приравнивая в каждом внутреннем узле x = x_i значения этих производных, вычисленных на концах рассматриваемого отрезка, получают (2n–2) уравнений

b_i+1 = b_i + 2h_ic_i + 3h d_i ; i=1,2,…,n–1; (35)

c_i+1 = c_i + 3h_id_i ; i=1,2,…,n–1 . (36)

Оставшиеся 2 уравнения получают из естественного предположения условия о нулевой кривизне этой функции на концах отрезка.

(37)

Система, составленная из (32) – (37), решается одним из методов решения СЛАУ.

Для упрощения машинных расчетов эта система уравнений приводится к более удобному виду посредством следующего алгоритма.

1. Из условия (32) можно сразу найти a_i.

2. Из (36) – (37) находят:

(38)

3. После подстановки (38) и (32) в (33) находят коэффициенты b_i.

b_i = ;

b_n = . (39)

4. Учитывая (38) и (39) из уравнения (35) исключаются d_i и b_i, тогда исходная система приводится к трехдиагональной матрице, содержащей только коэффициенты c_i. Получаем систему

h_i–1c_i–1 + 2(h_i–1 + h_i)c_i + h_ic_i+1 = 3( ), i=2,3,…,n. (40)

При этом c₁ = 0, c_n+1 = 0. Система (40) может быть решена методом прогонки. Зная c_i по (38) и (39), определяют b_i и d_i. Тогда кубический многочлен определяется для всех интервалов.

Пример составления системы (40). Пусть функция f(x) задана таблицей

i	0	1	2	3	4	5
x	0,1	0,15	0,19	0,25	0,28	0,30
y = f(x)	1,1052	1,1618	1,2092	1,2840	1,3231	0,3499
h		0,05	0,04	0,06	0,03	0,02

с₁ = 0;

0,05c₁+0,18c₂+0,04c₃ = 3

(коэффициент при c₂ получен следующим образом: 2(0,05+0,04) = 0,18;)

0,04c₂ + 0,2c₃ + 0,06c₄ = 3 ;

0,06c₃ + 0,18c₄ + 0,03c₅ = 3

0,03c₄ + 0,1c₅ = 3

c₆ = 0.

В результате получим систему относительно c₂  c₅:

 = .

Найдя c_i по (38), находят d_i и затем по (39) – b_i.