Інтерполювання із скінченими різницями за допомогою Excel

Означення 10. Якщо вузли інтерполяції х_k = х₀ + kh, де h – дійсна стала, k = 0,1,…,n, y_k = f (х_k), то таблицю

х						n
х0	y0	Δу0	Δ^2у0	Δ^3у0	…	Δ^nу0
х1	y1	Δу1	Δ^2у1	Δ^3у1	…
х2	y2	Δу2	Δ^2у2	Δ^3у2
	…	…	…
хn-2	yn-2	Δуn-2	Δ^2уn-2
хn-1	yn-1	Δуn-1
хn	yn

 

називають горизонтальною таблицею її скінчених різниць.

Задача 1. Знайти для функції f, заданої таблицею

 

i
xi		0,2	0,4	0,6	0,8
yi		0,1974	0,3805	0,5404	0,6747	0,7854

 

її скінчені різниці Δ²у₃ , Δ⁴у₃, Δ⁵у₀ .

Розв’язання. Побудуємо таблицю скінчених різниць, звідки й отримаємо шукані значення. Надамо чарункам електронної таблиці таких значень:

 

	A	B	C	D	E	F	G
			= B3 – B2	→	→	→	→
		0,1974	↓	↓	↓	↓
		0,3805	↓	↓	↓
		0,5404	↓	↓
		0,6747	↓
		0,7854

 

У рядку 1 указані порядки скінчених різниць у відповідному стовпці. У стовпці B тут різниці порядку 0, тобто значення функції f у вузлі інтерполяції, номер якого вказаний у стовпці А. Правила копіювання Excel дозволяють отримати у відповідних чарунках саме ті формули, які повинні там бути згідно з попереднім означенням: так у чарунці Е4 міститься Δ³у₂. В результаті отримаємо таку трикутну таблицю:

 

	A	B	C	D	E	F	G
			0,1974	-0,01428	-0,00891	0,006519	-0,0022
		0,1974	0,183111	-0,0232	-0,00239	0,004321
		0,3805	0,159913	-0,02559	0,001927
		0,5404	0,134321	-0,02366
		0,6747	0,110657
		0,7854

 

Тут Δ³у₂ ≈ 0,0019, Δ⁵у₀ ≈ – 0,0022, скінчена різниця Δ⁴у₃ не визначена умовами задачі.

Задача 2. Для функції f, заданої таблицею

 

i
xi		0,2	0,4	0,6	0,8		1,2
yi		0,1974	0,3805	0,5404	0,6747	0,7854	0,8761

 

знайти наближені значення функції в точках x₁* = 0,1 і x₂* = 1,1.

Розв’язання. Побудуємо спочатку таблицю скінчених різниць, за якою знайдемо інтерполяційний многочлен Ньютона. Оскільки дана в умові таблиця відрізняється від таблиці в умові задачі 1 тільки доданим шостим вузлом інтерполяції, то і таблицю скінчених різниць цієї задачі можна отримати з такої таблиці попередньої задачі, просто продовживши її копіювання:

	A	B	C	D	E	F	G	H
			= B3 – B2	→	→	→	→	→
		0,1974	↓	↓	↓	↓	↓
		0,3805	↓	↓	↓	↓
		0,5404	↓	↓	↓
		0,6747	↓	↓
		0,7854	↓
		0,8761

 

В результаті отримаємо таку трикутну таблицю:

 

	A	B	C	D	E	F	G	H
			0,1974	-0,01428	-0,00891	0,006519	-0,0022	-0,00038
		0,1974	0,183111	-0,0232	-0,00239	0,004321	-0,00258
		0,3805	0,159913	-0,02559	0,001927	0,001739
		0,5404	0,134321	-0,02366	0,003667
		0,6747	0,110657	-0,02
		0,7854	0,09066
		0,8761

 

Оскільки в даному випадку x₁* = 0,1 міститься на початку таблиці, тобто x₁* є [x₀; x₁], то тут доцільно буде скористатися формулою Ньютона дляінтерполюваннявперед, тобто першою інтерполяційною формулою Ньютона (3). У ній задіяні тільки скінчені різниці Δ^kу₀, які всі містяться у рядку 2 останньої таблиці. За умовою t = (x₁* – x₀)/h = (0,1 – 0)/0,2 = 0,5. Побудуємо таблицю для підрахунків L_k(x₁*) (k = 1,…,6). Надамо таких значень чарункам електронної таблиці:

 

	A	B	C	D	E	F	G	H
	k
	t-k	0,5	= 0,5 – C11	→	→	→	→	→
	t*…*(t-k+1)		= B13*B12	→	→	→	→	→
	k!		= B14*C11	→	→	→	→	→
	uk(x1*)	= B2*B13/B14	→	→	→	→	→	→
	Lk(x1*)	= СУММ($B$15:B15)	→	→	→	→	→	→

 

Тут у рядку 11 позначений номер k одночлена u_k(x) = Δ^kу₀ першого інтерполяційного многочлена Ньютона L_k(x). У чарунці В12 значення кроку t = 0,5 вузлів інтерполяції. У рядку 13 підраховуються множники t(t – 1)(t – 2)…(t – k + 1) з номером k, позначеним у рядку 11. У рядку 14 підраховуються значення факторіалу k! з відповідним для цього стовпця номером k. У рядку 15 знаходимо значення одночлена u_k(x) з відповідним стовпцю номером k у заданій точці x₁*. Нарешті у рядку 16 послідовно знаходимо суми одночленів u_k(x₁*), шукане значення L₆(x₁*) дістанемо у чарунці H16. В результаті отримаємо таку таблицю:

 

	A	B	C	D	E	F	G	H
	k
	t-k	0,5	-0,5	-1,5	-2,5	-3,5	-4,5	-5,5
	t*…*(t-k+1)		0,5	-0,25	0,375	-0,9375	3,28125	-14,7656
	k!
	uk(x1*)		0,098698	0,001786	-0,00056	-0,00025	-6E-05	7,89E-06
	Lk(x1*)		0,098698	0,100483	0,099926	0,099672	0,099612	0,099619

 

Отже, відповідь f (0,1) ≈ L₆(0,1) ≈ 0,099619. Щодо x₂* = 1,1, то це значення міститься ближче до кінця відрізкаінтерполювання, x₂* є [x₅; x₆]. Тож тут доцільно скористатися формулою Ньютона дляінтерполювання назад, тобто другою інтерполяційною формулою Ньютона (4). У ній задіяні тільки скінчені різниці Δ^kу_n-k (k = 1,…,n), які за означенням 9 всі містяться на діагоналі горизонтальної таблиці скінчених різниць. Оскільки зручніше для подальших підрахунків тримати всі ці значення на горизонталі, то таблицю скінчених різниць дещо переформуємо у порівнянні з означенням 10:

 

					n – 1	n
y0
y1	Δу0
y2	Δу1	Δ^2у0
…	…	…	…
yn-2	Δуn-3	Δ^2уn-4	Δ^3уn-5	…
yn-1	Δуn-2	Δ^2уn-3	Δ^3уn-4	…	Δ^n-1у0
yn	Δуn-1	Δ^2уn-2	Δ^3уn-3	…	Δ^n-1у1	Δ^nу0

 

Зауваження. Традиційно цю останню таблицю записують у такому форматі:

 

х	у	Δу	Δ^2у	Δ^3у	Δ^4у
х0	у0
		Δу0
х1	у1		Δ^2у0
		Δу1		Δ^3у0
х2	у2		Δ^2у1		Δ^4у0
		Δу3		Δ^3у1
х3	у3		Δ^2у2
		Δу4
х4	у4

Природно, що таку таблицю звичайно називають діагональною. Однак нам зручніше буде використовувати її для розрахунків у попередньому форматі.

Отже, надамо чарункам електронної таблиці таких значень:

 

	A	B	C	D	E	F	G	H
		0,1974	↑
		0,3805	↑	↑
		0,5404	↑	↑	↑
		0,6747	↑	↑	↑	↑
		0,7854	↑	↑	↑	↑	↑
		0,8761	= В25 – В24	→	→	→	→	→

 

В результаті отримаємо ту ж саму таблицю скінчених різниць у новому форматі:

 

	A	B	C	D	E	F	G	H
	0,2	0,1974	0,197396
	0,4	0,3805	0,183111	-0,01428
	0,6	0,5404	0,159913	-0,0232	-0,00891
	0,8	0,6747	0,134321	-0,02559	-0,00239	0,006519
		0,7854	0,110657	-0,02366	0,001927	0,004321	-0,0022
	1,2	0,8761	0,09066	-0,02	0,003667	0,001739	-0,00258	-0,00038

 

Побудуємо таблицю для підрахунку L_k(x₂*). Тут за умовою t = (x₂* – x₆)/h = (1,1 – 1,2)/0,2 = – 0,5. Надамо чарункам електронної таблиці таких значень:

 

	A	B	C	D	E	F	G	H
	k
	t + k	– 0,5	= – 0,5 + C27	→	→	→	→	→
	t*…*(t + k–1)		= B29*B28	→	→	→	→	→
	k!		= B30*C27	→	→	→	→	→
	uk(x2*)	= B25*B29/B30	→	→	→	→	→	→
	Lk(x2*)	= СУММ($B$32:B32)	→	→	→	→	→	→

Вона майже не відрізняється від копії таблиці для підрахунку L_k(x₁*) у нові чарунки: лише у чарунці В28 нове значення t = – 0,5, у формулі чарунки С28 заміняємо – на + і у чарунці В31 посилаємось на таблицю скінчених різниць у новому форматі. В результаті отримаємо таку таблицю:

 

	A	B	C	D	E	F	G	H
	k
	t+k	-0,5	0,5	1,5	2,5	3,5	4,5	5,5
	t*…(t+k-1)		-0,5	-0,25	-0,375	-0,9375	-3,28125	-14,7656
	k!
	uk(x2*)	0,876058	-0,04533	0,0025	-0,00023	-6,8E-05	7,06E-05	7,89E-06
	Lk(x2*)	0,876058	0,830728	0,833228	0,832999	0,832931	0,833001	0,833009

 

Шукане значення L₆(x₂*) у чарунці H32. Отже, відповідь f (1,1) ≈ L₆(1,1) ≈ 0,833.

Задача 3. Оцінити похибку інтерполяції у задачі 2, якщо вважати, що функція f (x) 9 разів неперервно диференційовна і задана таблицею

 

i
xi		0,2	0,4	0,6	0,8		1,2	1,4
yi		0,1974	0,3805	0,5404	0,6747	0,7854	0,8761	0,9505

Розв’язання. Згідно з теоремою 3 абсолютна похибка інтерполяції R₆(f,x) за першою інтерполяційною формулою Ньютона це R₆(f,x) = f (x) – L₆(x) ≈ t(t – 1)(t – 2)…(t – 7). Це останній одночлен u₇(x) першого інтерполяційного многочлена Ньютона L₇(x). Отже, для отримання похибки інтерполяції необхідно просто продовжити таблиці задачі 2 ще на восьмий вузол:

	A	B	C	D	E	F	G	H	I
			= B3 – B2	→	→	→	→	→	→
		0,1974	↓	↓	↓	↓	↓	↓
		0,3805	↓	↓	↓	↓	↓
		0,5404	↓	↓	↓	↓
		0,6747	↓	↓	↓
		0,7854	↓	↓
		0,8761	↓
		0,9595

 

В результаті отримаємо таку трикутну таблицю:

 

	A	B	C	D	E	F	G	H	I
			0,1974	-0,01428	-0,00891	0,006519	-0,0022	-0,00038	0,001386
		0,1974	0,183111	-0,0232	-0,00239	0,004321	-0,00258	0,001001
		0,3805	0,159913	-0,02559	0,001927	0,001739	0,00158
		0,5404	0,134321	-0,02366	0,003667	0,000159
		0,6747	0,110657	-0,02	0,003826
		0,7854	0,09066	-0,01617
		0,8761	0,074489
		0,9505

 

Побудуємо таблицю для підрахунку L₇(x₁*):

 

	A	B	C	D	E	F	G	H	I
	k
	t-k	0,5	= 0,5 – C11	→	→	→	→	→	→
	t*…*(t-k+1)		= B13*B12	→	→	→	→	→	→
	k!		= B14*C11	→	→	→	→	→	→
	uk(x1*)	= B2*B13/B14	→	→	→	→	→	→	→
	Lk(x1*)	= СУММ($B$15:B15)	→	→	→	→	→	→	→

 

В результаті отримаємо таку таблицю:

 

	A	B	C	D	E	F	G	H	I
	k
	t-k	0,5	-0,5	-1,5	-2,5	-3,5	-4,5	-5,5	-6,5
	t*…(t-k+1)		0,5	-0,25	0,375	-0,9375	3,28125	-14,7656	81,21094
	k!
	uk(x1*)		0,098698	0,001786	-0,00056	-0,00025	-6E-05	7,89E-06	2,23E-05
	Lk(x1*)		0,098698	0,100483	0,099926	0,099672	0,099612	0,099619	0,099642

 

Абсолютна похибка інтерполяції R₆(f,x) = u₇(x) міститься у чарунці I15: R₆(f,x) ≈ 2,23∙10^-5. Отже, результат отриманий з табличною точністю чотирьох вірних значущих цифр. Аналогічно абсолютна похибка інтерполяції R₆(f,x) за другою інтерполяційною формулою Ньютона це R₆(f,x) = f (x) – L₆(x) ≈ t(t +1)(t + 2)…(t + n) = u₇(x). Отже, для отримання похибки інтерполяції і тут треба продовжити відповідні таблиці на восьмий вузол:

 

	A	B	C	D	E	F	G	H	I
		0,1974	↑
		0,3805	↑	↑
		0,5404	↑	↑	↑
		0,6747	↑	↑	↑	↑
		0,7854	↑	↑	↑	↑	↑
		0,8761	↑	↑	↑	↑	↑	↑
		0,9595	= В26 – В25	→	→	→	→	→	→

 

В результаті отримаємо таку трикутну таблицю:

 

	A	B	C	D	E	F	G	H	I
	0,2	0,197396	0,197396
	0,4	0,380506	0,183111	-0,01428
	0,6	0,54042	0,159913	-0,0232	-0,00891
	0,8	0,674741	0,134321	-0,02559	-0,00239	0,006519
		0,785398	0,110657	-0,02366	0,001927	0,004321	-0,0022
	1,2	0,876058	0,09066	-0,02	0,003667	0,001739	-0,00258	-0,00038
	1,4	0,950547	0,074489	-0,01617	0,003826	0,000159	-0,00158	0,001002	0,001386

 

Побудуємо таблицю для підрахунку L₇(x₂*):

 

	A	B	C	D	E	F	G	H	I
	k
	t + k	– 0,5	= – 0,5 + C27	→	→	→	→	→	→
	t*…*(t + k–1)		= B29*B28	→	→	→	→	→	→
	k!		= B30*C27	→	→	→	→	→	→
	uk(x2*)	= B25*B29/B30	→	→	→	→	→	→	→
	Lk(x2*)	= СУММ($B$32:B32)	→	→	→	→	→	→	→

 

В результаті отримаємо:

 

	A	B	C	D	E	F	G	H	I
	k
	t+k	-0,5	0,5	1,5	2,5	3,5	4,5	5,5	6,5
	t*…(t+k-1)		-0,5	-0,25	-0,375	-0,9375	-3,28125	-14,7656	-81,2109
	k!
	uk(x2*)	0,950547	-0,03724	0,002021	-0,00024	-6,2E-06	4,32E-05	-2,1E-05	-2,2E-05
	Lk(x2*)	0,950547	0,913302	0,915324	0,915085	0,915078	0,915122	0,915101	0,915079

 

Абсолютна похибка інтерполяції R₆(f,x) = u₇(x) міститься у чарунці I31: R₆(f,x) ≈ 2,2∙10^-5. Отже, результат і тут отриманий з табличною точністю чотирьох вірних значущих цифр.

Означення 11. Емпіричною формулою називають наближену формулу, отриману інтерполяцієютаблиці експериментальних даних.

Задача 4. Побудувати емпіричну формулу для функції f, заданої таблицею

 

i
xi		0,1	0,2	0,3	0,4	0,5	0,6	0,7
yi		0,0997	0,1974	0,2915	0,3805	0,4636	0,5404	0,6107

Розв’язання. С початку складемо за даною таблицею електронну горизонтальну таблицю скінчених різниць:

 

	A	B	C	D	E	F	G
			0,099669	-0,00194	-0,00172	0,000378	7,13E-05
	0,1	0,0997	0,097727	-0,00367	-0,00135	0,000449	-1,4E-05
	0,2	0,1974	0,094061	-0,00501	-0,0009	0,000436	-7,1E-05
	0,3	0,2915	0,08905	-0,00591	-0,00046	0,000365
	0,4	0,3805	0,083141	-0,00637	-9,6E-05
	0,5	0,4636	0,076772	-0,00647
	0,6	0,5404	0,070306
	0,7	0,6107

 

Виявляється, що у стовпці G всі різниці п’ятого порядку з табличною точністю (чотири вірні значущі цифри) дорівнюють нулю. За означенням 7 скінчені різниці вищих порядків автоматично також нулі з такою точністю. Всі різниці четвертого порядку у стовпці F це стала 0,4 з точністю 10^-4. За наслідком з леми 2 §6 на відрізку інтерполяції f (x) ≈ L_m(x), де L_m(x) – інтерполяційний многочлен степеня m = 4 у будь – якій формі. Отже, тут у якості емпіричної формули природно застосувати першу інтерполяційну формулу Ньютона. В даному випадку за таблицею скінчених різниць це f (x) ≈ L₄(x) ≈ 0,1t – 0,0019 – 0,0017 + 0,0004 , де t = (х – x₀)/h = (х – 0)/0,1 = 10х, звідки L₄(x) ≈ х – 0,19 – 1,7 + 4 ≈ х – 0,1х(х – 0,1) – 0,3х(х – 0,1)(х – 0,2) + 0,173х(х – 0,1)(х – 0,2)(х – 0,3).

 

Інші методи інтерполювання

Для інтерполювання значення х функції, що знаходиться поблизу вузла х_i у середині таблиці з рівновіддаленими вузлами інтерполяції з кроком h, звичайно використовують формулу Стірлінга, якщо │t│= │x – х_i│/h ≤ 0,25 або формулу Бесселя, якщо 0,25 ≤│t│≤ 0,75. У цьому випадку зручно позначити цей найближчий вузол через х₀, а індекси інших вузлів рахувати так, як указано у наступній таблиці (де у_k = f (х_k)):

 

Індекс	Вузол	Функція
…	…	…	…	…	…	…	…	…
– 3	х-3	у-3	…	…	…	…	…	…
			Δу-3	…	…	…	…	…
– 2	х-2	у-2		Δ^2у-3	…	…	…	…
			Δу-2		Δ^3у-3	…	…	…
– 1	х-1	у-1		Δ^2у-2		Δ^4у-3	…	…
			Δу-1		Δ^3у-2		Δ^5у-3	…
	х0	у0		Δ^2у-1		Δ^4у-2		Δ^6у-3
			Δу0		Δ^3у-1		Δ^5у-2	…
	х1	у1		Δ^2у0		Δ^4у-1	…	…
			Δу1		Δ^3у0	…	…	…
	х2	у2		Δ^2у1	…	…	…	…
			Δу2	…	…	…	…	…
	х3	у3	…	…	…	…	…	…
…	…	…	…	…	…	…	…	…

Формула Стірлінга має вигляд L_n(x) = у₀ + t + Δ²у_-1 + + Δ⁴у_-2 + + Δ⁶у_-3 + … + Δ²ⁿу_-n, де t = (х – х₀)/h. Формула Бесселя має вигляд L_n(x) = + (t – )Δу₀ + + Δ³у_-1 + ∙ + Δ⁵у_-2 + + … + + + + Δ²ⁿ⁺¹у_-n, де t = (х – х₀)/h.

Якщо t = 1/2, то формула Бесселя значно спрощується і в такому разі називається формулою інтерполювання на середину: L_n(x) = – + ∙ – + … + (–1)ⁿ .

Зокрема при n = 1 отримаємо формулу квадратичної інтерполяції по Бесселю: L₁(x) = у₀ + tΔу₀ – (Δу₁ – Δу_-1). Електронні таблиці для підрахунків значень інтерполяційного многочлена аналогічні таким таблицям для формул Ньютона.

При інтерполюванні функцій з великою кількістю вузлів інтерполяційний многочлен має високій степінь, що спричиняє коливання на проміжках між вузлами. Щоб зменшити степінь інтерполяційного многочлена вузли інтерполювання можна розбити на групи і проводити інтерполювання по групах. Але в цьому разі порушуються аналітичні умови, з’являються точки розриву похідних. Позбутися цього недоліку можна за допомогою сплайнів – функцій, які гладко склеєні з різних частин многочленів заданого степеня.

Розіб’ємо заданий відрізок [a; b] на n частин точками а = х₀ < х₁ < … < х_n = b. Сплайном S_m(x) називається функція, що визначена на відрізку [a; b] і має на ньому неперервну похідну порядку m – 1, а на кожному відрізку [х_k; х_k+1] співпадає з деяким многочленом степеня не вище за m, причому хоча би на жодному з них степінь точно дорівнює m. Сплайн, який у вузлах х_k приймає ті ж самі значення, що й деяка функція f (x) називається інтерполяційним.

На практиці широко застосовують сплайни третього степеня (кубічні сплайни S₃(x)). Для побудови інтерполяційного кубічного сплайну розіб’ємо відрізок [a; b] на n рівних частин довжини h = (b – а)/ n. В цьому разі на відрізку [х_k; х_k+1] (k = 0,1,…n – 1) кубічний сплайн запишеться так:

S₃(x) = +

+ , (5)

де m_k, m_k+1 – деякі числа. Виявляється, що саме за такої формули S₃(х_k) = f (х_k), S₃(х_k+1) = f (х_k+1), S´₃(х_k) = m_k, S´₃(х_k+1) = m_k+1, що і визначає цю формулу. Отже, кубічний сплайн S₃(x) має неперервну першу похідну всюди на [a; b]. Тепер оберемо числа m_k так, щоби S₃(x) мав і другу неперервну похідну на [a; b]. Ця умова зводиться до умови рівності другої похідної зліва і справа у кожному вузлі х_k (k = 1,…,n – 1), вона приймає вигляд:

m_k-1 + 4m_k + m_k+1 = (f (х_k+1) – f (х_k)), (k = 1,…,n – 1). Це система n – 1 лінійних рівнянь відносно n + 1 невідомих m_k (k = 0,1,…,n). Для однозначного визначення невідомих додають ще дві умови: це умови обмежень на значення сплайна і його похідних на кінцях відрізка [a; b], які і називають крайовими. Існує декілька видів крайових умов, одна з найпоширеніших це S´₃(а) = f ´(а), S´₃(b) = f ´(b). Тоді система n + 1 лінійних рівнянь відносно n + 1 невідомих m_k (k = 0,1,…,n) має вигляд:

Матриця цієї системи не вироджена і тому вона має єдиний розв’язок. Розв’язати її можна методом Гаусса або будь – яким іншим методом розділу 3. Знайшовши m_k, побудувати інтерполяційний кубічний сплайн можна за формулою (5) на кожному відрізку [х_k; х_k+1] (k = 0,1,…,n) окремо.

Доведено, що кубічна сплайн – функція – це єдина функція Р(х) з мінімальною кривиною серед усіх функцій, які інтерполюють дану функцію f (x), тобто у всіх вузлах інтерполяції задовольняють умови Р(x₀) = f (x₀), …, Р(x_n) = f (x_n), і мають квадратично інтегровну другу похідну. Отже, в цьому розумінні кубічний сплайн найкраща з функцій, які інтерполюють дану функцію.

 

 

Питання, тести

 

1. Задача інтерполювання функції з n вузлами інтерполяції має єдиний розв’язок, якщо інтерполююча функція Р(х) – це

 

А	кубічний сплайн
Б	многочлен степеня n – 1
В	многочлен степеня n
Г	ряд Тейлора даної функції

 

1. Інтерполяційний многочлен Лагранжа L_n(x) вузлів інтерполяції має

 

А	Б	В	Г
n	n + 1	n – 1	будь – яку кількість

 

1. Інтерполяційна формула Лагранжа f (x) ≈ L₂(x) має вигляд