Uwaga. Jeżeli wprowadzimy oznaczenie yi = f (x0 + ih) , dla i = 1,2,..., n , to
k
k
j
"k yi = ìÅ‚ ÷Å‚ .
"(-1) Å" ëÅ‚ öÅ‚ Å" yi+k - j
ìÅ‚ ÷Å‚
j
j=0
íÅ‚ Å‚Å‚
Wzory interpolacyjne Newtona  pierwszy wzór interpolacyjny.
Załóżmy, że "h > 0 "n " N xi - xi-1 = h. Wówczas "i " N xi = x0 + i Å" h. Przyjmijmy oznaczenie
x - x0
q = . Wez
my pod uwagę następujący układ funkcji bazowych
h
Åš0 (x) = 1, Åš1(x) = q, Åš2 (x) = q(q -1),...,Åšn (x) = q(q -1)(q - 2)....(q - n +1).
Wówczas wielomian interpolacyjny w można zapisać w postaci:
n
" x " R w(x) = Å" Åši (x) = b0 + b1q + b2q(q -1) +...+ bnq(q -1)...(q - n +1).
"bi
i=0
x0 - x0
Ponieważ dla x = x0 q = = 0,
h
x1 - x0 h
x = x1 q = = = 1,
h h
x2 - x0 2h
x = x2 q = = = 2,
h h
& & & & & & & & & & & & .
xn - x0 nh
x = xn q = = = n,
h h
więc współczynniki wielomianu interpolacyjnego w można wyznaczyć rozwiązując układ równań:
w(x0 ) = b0 = y0
Å„Å‚
ôÅ‚w(x ) = b0 + b1
= y1
1
ôÅ‚
ôÅ‚
= y2
òÅ‚w(x ) = b0 + 2b1 + 2b2
2
ôÅ‚
ôÅ‚.......................................................................................................
ôÅ‚w(xn ) = b0 + na1 + n(n -1)b2 + n(n -1)(n - 2)b3 +....+ n!bn = yn
ół
Z układu tego wynika, że
1 1 1
b0 = y0 , b1 = y1 - b0 = y1 - y0 = "y0, b2 = (y2 - 2b1 - b0 ) = [(y2 - b1 - b0 ) - b1] = ["y1 - "y0 ] =
2 2 2
"2 y0
= , itd.
2!
"i y0
Zatem "i " P(n) bn = . Stąd wynika, że pierwszy wzór interpolacyjny Newtona ma postać:
n!
y0 "y0 "2 y0 "n y0
w(x) = + q + q(q -1) +...+ q(q -1)...(q - n +1).
0! 1! 2! n!
Algorytm obliczania wartości funkcji za pomocą wielomianu interpolacyjnego Newtona
Dane: n - całkowite, h, x0, q  rzeczywiste, y[n]  wektor
Podaj wartości n, x0, h, x
od i = 0 do n podaj wartości yi0
od i = 1 do n
od j = 0 do n-i wykonuj
y = y - y tablica różnic skończonych
ji j+1,i-1 j,i-1
x - x0
q =
h
od i = 0 do n wykonuj
bi = y0i
od j = 0 do n wykonuj
bi
bi = (q - j)
j +1
w = 0
od i = 0 do n wykonuj
w = w + bi
Drukuj w.
Numeryczne całkowanie funkcji.
b
Numeryczne obliczanie wartości całki f (x)dx polega na poszukiwaniu jej wartości w postaci
+"
a
kombinacji liniowej wartości funkcji podcałkowej f(x) i jej pochodnych do pewnego ustalonego rzędu
włącznie w ustalonych punktach węzłowych.
m
Niech n, m " N, f :[a,b] ‚" R R bÄ™dzie funkcjÄ… klasy C ([a,b]) , zaÅ›
"i " P(m) "k " P(n) xik "[a,b] ustalonymi punktami. Szukamy takich liczb Aik dla i =1,2,& ,m;
b
n m
(i-1)
k =0,1,& ,n, aby f (x)dx H" Q( f ) = Aik Å" f (xik ).
" "
+"
k=0 i=1
a
Def. Powyższe przedstawienie wartości całki nazywamy kwadraturą liniową, liczby Aik nazywamy
współczynnikami kwadratury zaś punkty xik nazywamy węzłami kwadratury.
Rozważamy zadanie polegające na takim doborze wartości współczynników i węzłów, aby
kwadratura była dokładna w przypadku, gdy funkcją f jest wielomian określonego stopnia.
Rozważymy obecnie przypadek, gdy znane są jedynie wartości funkcji f w węzłach
x0 , x1,..., xn. Załóżmy, że węzły są uporządkowane rosnąco i są parami różne, tzn. x0 < x1 <...< xn.
n
Rozważana kwadratura ma zatem postać: Q( f ) = Ak Å" f (xk ).
"
k=0
Def. Mówimy, że kwadratura Q jest rzędu r " N , jeśli dla dowolnego wielomianu w stopnia
b b
mniejszego od r : Q(w) = w(x)dx oraz istnieje wielomian w stopnia r taki, że Q(w) `" w(x)dx.
+" +"
a a
Kwadratury Newtona  Cotesa.
b
Def. Kwadraturami Newtona  Cotesa przybliżającymi wartość całki f (x)dx nazywamy kwadratury
+"
a
otrzymane poprzez zastąpienie funkcji podcałkowej jej wielomianem interpolacyjnym Lagrange a
b - a
opartym na równoodlegÅ‚ych wÄ™zÅ‚ach "i " P(n) xi = a + i Å" h, h = .
n
b
W tym przypadku Q( f ) = Ln (x)dx.
+"
a
Przypomnijmy, że po wprowadzeniu oznaczenia x = a + t Å" h wielomian Ln można zapisać w postaci:
n
b - a
Ln (a + th) = f (xi ) Å" , gdzie "i " P(n) xi = a + i Å" h, h = .
" "t - j
i - j n
i=0 j=0
j`"i
Zatem w tym przypadku
b b n
n n n n n n n
x - x
j
Q( f ) = f (xi )Å" f (xi )Å" Ai Å" f (xi ),
" " " "t - jd(a + th) =" "t - j Å" h Å" dt ="
+"i=0 j=0 xi - x dx = f (xi )Å"+" +"
i - j i - j
i=0 j=0 i=0 j=0 i=0
j
a a 0
j`"i j`"i j`"i
n
n
gdzie (*) Ai = h Å"
"t - jdt.
+"
i - j
j=0
0
j`"i
Lemat. "i " P(n) Ai = An-i.
Przykłady najprostszych kwadratur Newtona  Cotesa.
b
1. Dla n = 0 f (x)dx H" Q( f ) = h Å" f (a) . Ten sposób obliczania caÅ‚ki nosi nazwÄ™ metody
+"
a
prostokątów.
Twierdzenie. Jeżeli f jest funkcją klasy C1([a,b]) , to błąd metody prostokątów wynosi
(b - a)2
2
R( f ) = Å" f (¾ ) , gdzie ¾ " (a,b).
2
Wzór ten możemy zastosować w każdym z przedziałów [xi , xi+1] dla i = 0,1,& ,n-1.
n-1
Otrzymamy wówczas, że Q( f ) = f (xi ) Å" h.
"
i=0
2. Dla n = 1 węzłami kwadratury są krańce przedziału całkowania. Ze wzoru (*) obliczamy
1
2
t -1 t b - a
współczynniki kwadratury A0 = (b - a) Å" dt = (b - a) Å"[- + t]1 = ,
0
+"
0 -1 2 2
0
1
2
t - 0 t b - a
A1 = (b - a) Å" dt = (b - a) Å"[ ]1 = . Kwadratura N  C dla n = 1 ma wiÄ™c postać:
0
+"1- 0
2 2
0
1
b - a
Q( f ) = Ai Å" f (xi ) = Å" ( f (a) + f (b)). Ten sposób obliczania caÅ‚ki nosi nazwÄ™ metody
"
2
i=0
trapezów.
2
Twierdzenie. Jeżeli f " C ([a,b]), to błąd metody trapezów wyraża się wzorem
(b - a)3
2 2
R( f ) = Å" f (¾ ) , gdzie ¾ " (a,b).
12
Wzór ten możemy zastosować w każdym z przedziałów [xi , xi+1] dla i = 0,1,& ,n-1.
n-1
xi+1 - xi y0 + yn
ëÅ‚ öÅ‚
Otrzymamy wówczas, że Q( f ) = Å"( f (xi-1) + f (xi )) = h Å" + y1 +...+ yn-1 .
ìÅ‚ ÷Å‚
"
2 2
i=0 íÅ‚ Å‚Å‚
a + b
3. Dla n = 2 węzłami kwadratury są x0 = a, x1 = , x2 = b. Ze wzoru (*) wynika, że
2
współczynniki kwadratury w tym przypadku są równe odpowiednio
2
2
3 2
îÅ‚ Å‚Å‚ - a 8 12
b - a (t -1)(t - 2) b - a t 3t b b - a
îÅ‚
A0 = Å" dt = Å" - + 2tśł = Å" - + 4Å‚Å‚ = ,
ïÅ‚
+" ïÅ‚3 śł
2 (0 -1)(0 - 2) 4 3 2
ðÅ‚ ûÅ‚0 4 ðÅ‚ 2 ûÅ‚ 6
0
2
b - a (t - 0)(t - 2) 4(b - a)
A1 = Å" dt = , A2 = A0.
+"
2 (1- 0)(1- 2) 6
0
Otrzymana w tym przypadku kwadratura ma postać
b - a a + b b - a
îÅ‚
Q( f ) = f (a) + 4 f ( ) + f (b)Å‚Å‚ = Å" (y0 + 4y1 + y2 ) i nazywana jest wzorem
ïÅ‚ śł
6 2 6
ðÅ‚ ûÅ‚
parabol lub wzorem Simpsona.
4
Twierdzenie. Rząd metody Simpsona wynosi 4. Ponadto, jeżeli f " C ([a,b]), to błąd metody
5
1 b - a
ëÅ‚ öÅ‚
(4)
Simpsona wyraża siÄ™ wzorem R( f ) = Å" f (¾ ) , gdzie ¾ " (a,b).
ìÅ‚ ÷Å‚
90 2
íÅ‚ Å‚Å‚
Wzór ten możemy zastosować w każdym z przedziałów [xi , xi+2 ] dla i = 0,1,& ,2k.
(oczywiście, gdy n = 2k ). Otrzymamy wówczas, że
b - a b - a
Q( f ) = Å"(y0 + 4y1 + 2y2 + 4y3 +...+ 2yn-2 + 4yn-1 + yn ). ( bo, = h )
3n 2n
Błędy poszczególnych metod dla węzłów x0 , x1,..., xn wyrażają się wzorami:
(b - a)2
2
- dla metody prostokÄ…tów R( f ) = Å" f (¾ ) , gdzie ¾ " (a,b).
2n
(b - a)3
2 2
- dla metody trapezów R( f ) = Å" f (¾ ) , gdzie ¾ " (a,b).
12n2
5
1 b - a
ëÅ‚ öÅ‚
(4)
- dla metody Simpsona R( f ) = Å" f (¾ ) , gdzie ¾ " (a,b).
ìÅ‚ ÷Å‚
180n4 2
íÅ‚ Å‚Å‚
Twierdzenie. Kwadratury Newtona  Cotesa oparte na n+1 węzłach są rzędu n+2 dla n  parzystych,
zaś rzędu n+1 dla n  nieparzystych.
Algorytmy dla poszczególnych metod:
Metoda prostokątów Metoda trapezów Metoda Simpsona
b - a b - a f (a) + f (b) b - a
h = , c = 0 h = , c = h = , c = f (a) + f (b), p = 1
n n 2 n
Dla i = 0, 1, & , n-1 wykonuj Dla i = 1, 2, & , n-1 wykonuj Dla i = 1, 2, & , n-1 wykonuj
x=a+ih x=a+ih x=a+ih
c=c+f(x) c=c+f(x) jeżeli p = 1 to p:=0, c:=c+4f(x)
Q(f)=hc Q(f)=hc else p:=1, c:=c+2f(x)
hc
Q(f)=
3
Numeryczne różniczkowanie funkcji
Niech f:[a,b]R będzie funkcją różniczkowalną na przedziale (a,b), x0 , x1,..., xn "[a,b], gdzie
x0 < x1 <...< xn , zaś y0 = f (x0 ), y1 = f (x1),..., yn = f (xn ) będą znanymi wartościami funkcji f.
Szukamy wartości i - tej pochodnej funkcji f w punkcie x "[x0 , xn ]. Rozwiązanie tego problemu
polega na interpolacji funkcji f na przedziale [a,b] wielomianem interpolacyjnym w, a następnie
przyjęciu, że wartość i  tej pochodnej funkcji f w punkcie x równa jest w przybliżeniu wartości i  tej
i i
d f d w
pochodnej wielomianu w tym punkcie (x) H" (x) dla x "[x0 , xn ].
dxi dxi
Różniczkowanie funkcji w oparciu o wielomian interpolacyjny Newtona
Załóżmy, że wÄ™zÅ‚y x0 , x1,..., xn sÄ… równoodlegÅ‚e tzn. " h > 0 "k " P(n) xk = x0 + k Å" h. Oznaczmy
x - x0
przez q = i rozważmy pierwszy wzór interpolacyjny Newtona
h
y0 "y0 "2 y0 "n y0
w(x) = + q + q(q -1) +...+ q(q -1)...(q - n +1).
0! 1! 2! n!
dw dw dq dq 1
KorzystajÄ…c z reguÅ‚ różniczkowania funkcji zÅ‚ożonej otrzymamy, że = Å" . Ponieważ =
dx dq dx dx h
dw 2q -1 3q2 - 6q + 2
i = "y0 + "2 y0 Å" + "3 y0 Å" +... , wiÄ™c
dq 2! 3!
ëÅ‚ öÅ‚
df dw 1 2q -1 3q2 - 6q + 2
ìÅ‚ ÷Å‚
" x "[x0 , xn ] (x) H" (x) = Å" + "2 y0 Å" + "3 y0 Å" +...÷Å‚.
ìÅ‚"y0
dx dx h 2! 3!
íÅ‚ Å‚Å‚
Stąd w szczególności otrzymamy, że
ëÅ‚ öÅ‚
df dw 1 "2 y0 2"3 y0
ìÅ‚ ÷Å‚
- dla q = 0 (x0 ) H" (x0 ) = Å" - +
ìÅ‚"y0 2! 3! +...÷Å‚,
dx dx h
íÅ‚ Å‚Å‚
ëÅ‚ öÅ‚
df dw 1 "2 y0 "3 y0
ìÅ‚
- dla q = 1 (x1) H" (x1) = Å" + - +...÷Å‚,
ìÅ‚"y0 2! 3! ÷Å‚
dx dx h
íÅ‚ Å‚Å‚
ëÅ‚ öÅ‚
df dw 1 3"2 y0 2"3 y0
ìÅ‚ ÷Å‚
- dla q = 2 (x2 ) H" (x2 ) = Å" +
ìÅ‚"y0 2! + 3! +...÷Å‚.
dx dx h
íÅ‚ Å‚Å‚
Powyższe wzory określają wartości pochodnych funkcji f w punktach węzłowych.
Rozważmy przypadek, gdy znane są wartości funkcji jedynie w trzech węzłach równoodległych
y0 = f (x0 ), y1 = f (x0 + h), y2 = f (x0 + 2h). StosujÄ…c oznaczenie x1 = x0 + h, x2 = x0 + 2h
stwierdzamy, że możemy obliczyć jedynie dwie różnice "y0 = y1 - y0 , oraz
"2 y0 = (y2 - y1) - (y1 - y0 ) = y2 - 2y1 + y0. W tym przypadku naturalnym jest pominięcie wyrazów
zawierających różnice rzędu wyższego niż 2. Otrzymamy zatem, że
öÅ‚
df dw 1 ëÅ‚ 2q -1
ìÅ‚
" x "[x0 , xn ] (x) H" (x) = Å" y1 - y0 + ( y2 - 2y1 + y0 ) Å" ÷Å‚
ìÅ‚ ÷Å‚,
dx dx h 2!
íÅ‚ Å‚Å‚
y2 - 2y1 + y0
df dw 1 ëÅ‚ öÅ‚ 1
(x0 ) H" (x0 ) = Å" y1 - y0 - = Å" (-3y0 + 4y1 - y2 ),
ìÅ‚ ÷Å‚
dx dx h 2 2h
íÅ‚ Å‚Å‚
y2 - 2y1 + y0
df dw 1 ëÅ‚ öÅ‚ 1
(x1) H" (x1) = Å" y1 - y0 + = Å" ( y2 - y0 ),
ìÅ‚ ÷Å‚
dx dx h 2 2h
íÅ‚ Å‚Å‚
3( y2 - 2y1 + y0
df dw 1 ëÅ‚ öÅ‚ 1
(x2 ) H" (x2 ) = Å" y1 - y0 + = Å" (y0 - 4y1 + 3y2 ).
ìÅ‚ ÷Å‚
dx dx h 2 2h
íÅ‚ Å‚Å‚
Wyznaczanie pochodnej rzędu drugiego
Drugiego rzędu pochodna wielomianu interpolacyjnego Newtona ma postać:
2 2
2 2 2 2 2 2
d w d dw dq d w dq dw d q d w dq 1 d w dq 1 d q
ëÅ‚ öÅ‚ ëÅ‚ öÅ‚
= ( Å" ) = Å" + Å" = Å" = Å" , bo = i = 0. StÄ…d
ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚
dx2 dx dq dx dq2 dx dq dx2 dq2 dx h2 dq2 dx h dx2
íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚
2
ëÅ‚ öÅ‚
d w 1 2"2 y0 - 6 1
6q
2 2 ìÅ‚ ÷Å‚
otrzymamy, że f (x) H" (x) = Å" + Å" "3 y0 +...÷Å‚ = Å"("2 y0 + (q -1) Å" "3 y0 +...).
dx2 h2 ìÅ‚ 2! 3! h2
íÅ‚ Å‚Å‚
Wstawiając kolejno wartości q = 0, q = 1, q = 2, otrzymamy wartości pochodnej rzędu drugiego w
1 1
2 2 2 2 2 2 2 2
wÄ™zÅ‚ach f (x0 ) H" w (x0 ) = Å" ("2 y0 - "3 y0 +...), f (x1) H" w (x1) = Å" ("2 y0 +...),
h2 h2
1
2 2 2 2
f (x2 ) H" w (x2 ) = Å" ("2 y0 + "3 y0 +...).
h2
Przy znajomości wartości funkcji tylko w trzech węzłach, istnieje możliwość obliczenia
1
2 2
pochodnej rzędu drugiego funkcji f jedynie w węz
le x1: f (x1) H" Å"(y2 - 2y1 + y0 ).
h2
W pozostałych węzłach otrzymalibyśmy ten sam wynik, który niewiele ma wspólnego z
rzeczywistością.
Znając wartości funkcji w czterech węzłach, możemy sensownie obliczać wartości pochodnej
we wszystkich węzłach. I tak np:
1
2 2
f (x0 ) H" Å"(2y0 - 5y1 + 4y2 - y3),
h2
1
2 2
f (x1) H" Å"(- 2y1 + 2y2 ),
h2
1
2 2
f (x2 ) H" Å"(y1 - 2y0 + y3).
h2
Uwaga: Można pokazać, że błędy metod interpolacyjnych są proporcjonalne do wartości kroku h,
natomiast błędy zaokrągleń są odwrotnie proporcjonalne do wartości hk , gdzie k jest maksymalnym
wykładnikiem we wzorze interpolacyjnym. Zatem zmniejszanie kroku h powoduje z jednej strony
zmniejszanie błędu interpolacji, a z drugiej strony błędy zaokrąglenia powiększają się. Zatem w
przypadku, gdy mamy możliwość liczyć wartości funkcji w dowolnym punkcie, nie należy przesadzać
ze zmniejszaniem wartości kroku h.