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



Мы поможем в написании ваших работ!


Мы поможем в написании ваших работ!



Мы поможем в написании ваших работ!


ЗНАЕТЕ ЛИ ВЫ?

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



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

х   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

 

її скінчені різниці Δ2у3 , Δ4у3 , Δ5у0 .

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

 

  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 міститься Δ3у2. В результаті отримаємо таку трикутну таблицю:

 

  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          

 

Тут Δ3у2 ≈ 0,0019, Δ5у0 ≈ – 0,0022, скінчена різниця Δ4у3 не визначена умовами задачі.

Задача 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

 

знайти наближені значення функції в точках x1* = 0,1 і x2* = 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            

 

Оскільки в даному випадку x1* = 0,1 міститься на початку таблиці, тобто x1* є [x0 ; x1], то тут доцільно буде скористатися формулою Ньютона дляінтерполюваннявперед, тобто першою інтерполяційною формулою Ньютона (3). У ній задіяні тільки скінчені різниці Δkу0 , які всі містяться у рядку 2 останньої таблиці. За умовою t = (x1* – x0)/h = (0,1 – 0)/0,2 = 0,5. Побудуємо таблицю для підрахунків Lk(x1*) (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 одночлена uk(x) = Δkу0 першого інтерполяційного многочлена Ньютона Lk(x). У чарунці В12 значення кроку t = 0,5 вузлів інтерполяції. У рядку 13 підраховуються множники t(t – 1)(t – 2)…(t – k + 1) з номером k, позначеним у рядку 11. У рядку 14 підраховуються значення факторіалу k! з відповідним для цього стовпця номером k. У рядку 15 знаходимо значення одночлена uk(x) з відповідним стовпцю номером k у заданій точці x1*. Нарешті у рядку 16 послідовно знаходимо суми одночленів uk(x1*), шукане значення L6(x1*) дістанемо у чарунці 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) ≈ L6(0,1) ≈ 0,099619. Щодо x2* = 1,1 , то це значення міститься ближче до кінця відрізкаінтерполювання, x2* є [x5 ; x6]. Тож тут доцільно скористатися формулою Ньютона дляінтерполювання назад, тобто другою інтерполяційною формулою Ньютона (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

 

Побудуємо таблицю для підрахунку Lk(x2*). Тут за умовою t = (x2* – x6)/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)

Вона майже не відрізняється від копії таблиці для підрахунку Lk(x1*) у нові чарунки: лише у чарунці В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

 

Шукане значення L6(x2*) у чарунці H32. Отже, відповідь f(1,1) ≈ L6(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 абсолютна похибка інтерполяції R6(f,x) за першою інтерполяційною формулою Ньютона це R6(f,x) = f(x) – L6(x) ≈ t(t – 1)(t – 2)…(t – 7). Це останній одночлен u7(x) першого інтерполяційного многочлена Ньютона L7(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              

 

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

 

  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

 

Абсолютна похибка інтерполяції R6(f,x) = u7(x) міститься у чарунці I15: R6(f,x) ≈ 2,23∙10-5. Отже, результат отриманий з табличною точністю чотирьох вірних значущих цифр. Аналогічно абсолютна похибка інтерполяції R6(f,x) за другою інтерполяційною формулою Ньютона це R6(f,x) = f(x) – L6(x) ≈ t(t +1)(t + 2)…(t + n) = u7(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

 

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

 

  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

 

Абсолютна похибка інтерполяції R6(f,x) = u7(x) міститься у чарунці I31: R6(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) ≈ Lm(x), де Lm(x) – інтерполяційний многочлен степеня m = 4 у будь – якій формі. Отже, тут у якості емпіричної формули природно застосувати першу інтерполяційну формулу Ньютона. В даному випадку за таблицею скінчених різниць це f(x) ≈ L4(x) ≈ 0,1t – 0,0019 – 0,0017 + 0,0004 , де t = (х – x0)/h = (х – 0)/0,1 = 10х, звідки L4(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. У цьому випадку зручно позначити цей найближчий вузол через х0 , а індекси інших вузлів рахувати так, як указано у наступній таблиці ( де уk = fk) ) :

 

Індекс Вузол Функція            
– 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

Формула Стірлінга має вигляд Ln(x) = у0 + t + Δ2у-1 + + Δ4у-2 + + Δ6у-3 + … + Δ2nу-n , де t = (х – х0)/h. Формула Бесселя має вигляд Ln(x) = + (t – )Δу0 + + Δ3у-1 + + Δ5у-2 + + … + + + + Δ2n+1у-n , де t = (х – х0)/h.

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

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

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

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

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

S3(x) = +

+ , (5)

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

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

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

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

 

 

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

 

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

 

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

 

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

 

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

 

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

 



Последнее изменение этой страницы: 2016-04-23; Нарушение авторского права страницы; Мы поможем в написании вашей работы!

infopedia.su Все материалы представленные на сайте исключительно с целью ознакомления читателями и не преследуют коммерческих целей или нарушение авторских прав. Обратная связь - 35.170.64.36 (0.018 с.)