7 - kwadratura

1

Kwadratury

1.1

Sformulowanie zadania

Naszym celem jest przyblizone obliczanie calek. Rozpoczniemy od przedstawienia metod numerycznego wyznaczania calek oznaczonych.

De nicja

Calka oznaczona Riemanna funkcji f na [a; b] nazywamy:

Z

b a

f (x)dx = nlim !1 Sn

(1)

gdzie: f - funkcja P okreslona i ograniczona na odcinku [a; b]; Sn = ni=0 (xi+1 xi )f (^xi ), dla dowolnego podzialu odcinka [a; b] ( a = x0 < x1 < ::: < xn+1 = b ) takiego, ze max0in j xi+1 xi j! 0 przy n ! 1 i dowolnych punktow posrednich x^i 2 [xi ; xi+1 ]. Funkcjonal liniowy, ktorym jest calka, oznaczamy przez I:

I (f ) =

Z

b

a

f (x)dx :

(2)

Przyblizone obliczanie calek mozna traktowac jako aproksymacje funkcjonalu I jakimis prostrzymi do obliczania funkcjonalami. W rachunku numerycznym musimy miec mozliwosc obliczania ich wartosci za pomoca skonczonej liczby dzialan arytmetycznych. Funkcjonalem tego typu jest funkcjonaly Q, tzw. kwadratury liniowe, postaci:

Q(f ) =

n0 X i=0

Ai;0 f (xi;0 ) +

n1 X i=0

X Ai;1 f 0(xi;1 ) + ::: + Ai;k f (k) (xi;k ) nk

i=0

(3)

gdzie: Ai;j - wspolczynniki kwadratury Q, xi;j - wezly kwadratury Q. Bedziemy zajmowali sie calkami postaci:

Ip (f ) =

Z

b a

p(x)f (x)dx :

(4)

gdzie funcja wagowa p jest nieujemna na odcinku [a; b], zeruje sie w nim w skonczonej liczbie punktow i Rb jest calkowalna, tzn. a p(x)dx < 1. Najczesciej stosuje sie kwadratury korzystajace jedynie z wartosci funkcji f :

Q(f ) = De nicja Reszte kwadratury

de niujemy jako:

n X

Ai f (xi ) :

(5)

R(f ) = Ip (f ) Q(f )

(6)

i=0

De nicja

Mowimy, ze kwadratura Q jest rzedu n jesli jest dokladna dla wszystkich wielomianow stopnia mniejszego od n, Q(w) = Ip (w) dla w 2 Wn 1 , oraz istnieje wielomian wn stopnia n, dla ktorego Q(wn ) 6= Ip (wn ).

1

1.2

Kwadratury interpolacyjne

Jedna z metod otrzymywania kwadratur postaci (5) jest calkowanie wielomianu lub funkcji sklejanej interpolujacej funkcje podcalkowa. Dokladniej poprzez kwadrautury interpolacyjne rozumie sie kwadratury otrzymane przez calkowanie wielomianow interpolacyjnych Hermite'a (w szczegolnosci Lagrange'a) funkcji podcalkowej f.

Lemat

Kwadratury interpolacyjne oparte na wezlach o lacznej krotnosci n + 1 sa rzedu conajmniej n + 1. Najprostszym przykladem kwadratury interpolacyjnej jest przedstawiony ponizej wzor trapezow.

1.2.1 Wzor trapezow Funkcje podcalkowa f przyblizamy wielomianem interpolacyjnym Lagrange'a opartym na wezlach a i b. f (x)  L1(x) = f (a)

x b x a + f (b) a b b a

(7)

Calkujac wielomian L1 (x) otrzymujemy kwadrature postaci:

Q(f ) = I (L1 ) =



Z b

x a b a b a x b f (a) + f (b) dx = f (a) + f (b) a b b a 2 2

a

(8)

Rb

Dla funkcji f dodatniej na [a; b], przyblizanie calki a f (x)dx powyzsza kwadratura mozna interpretowac geometrycznie jako przblizenie pola trapezu krzywoliniowego polem trapezu o podstawach f (a) i f (b). Dlatego wzor (8) nosi nazwe wzoru trapezow.

1.2.2 Funkcja sklejana Funkcje podcalkowa f mozna tex przyblizyc funkcja sklejana interpolujaca S1 , stopnia pierwszego, z wezlami a = x0 < x1 ; ::: < xN = b. W kazdym z przedzialow [xi ; xi+1 ]; i = 0; 1; :::; N 1, funkcja S1 jest okreslona wzorem

S1 (x) = f (xi ) +

f (xi+1 f (xi ) (x xi ) xi+1 xi

Calkujac S1 dostajemy kwadrature postaci

Q(f ) = I (S1 ) =

NX1 i=0

Z

xi+1 xi

S1 (x)dx =

NX1 i=0

x

i+1

Rb

2

xi

f (x i ) +

(9)

xi+1

2

xi



f (xi+1 )

(10)

Dla funkcji f dodatniej na [a; b], przyblizanie calki a f (x)dx powyzsza kwadratura mozna interpretowac geometrycznie jako przblizenie pola trapezu krzywoliniowego suma pol trapezow o wysokosciach xi+1 xi i podstawach f (xi ) i f (xi+1 ).

UWAGA

Kwadratura korzystajaca z funkcji sklejanych nie jest kwadratura interpolacyjna dla N > 1! Nie jest dla niej prawdziwy Lemat 1.

2

1.3

Kwadratury Newtona-Cotesa

Kwadraturami Newtona-Cotesa

R przyblizajacymi ab f (x) sa nazywane kwadratury Q(f ) = I (Ln )

(11)

gdzie Ln jest wielomianem interpolacyjnym Lagrange'a funkcji f opartym na rownoodleglych wezlach

x0 = a; x1 = a + x; : : : ; xn = a + nh = b

(12)

Stosujac podstawienie x = a + th mozemy zapisac wielomian Ln w postaci n X

Ln (x) = Ln(a + th) =

i=0

f (xi )

n Y t j j =0;j 6=i

gdzie

xi = a + ih; i = 0; 1; : : : ; n; h = Calkujac prawa strone (13) otrzymujemy kwadrature

Q(f ) = ze wspolczynnikami Ai okreslonymi wzorem

Ai = h

Z

n X

b a n

(14)

Ai f (xi )

(15)

n Y t j dt i j

(16)

i=0

n

(13)

i j

j =0;j 6=i

0

latwo sprawdzic, ze Ai = An i .

1.3.1 n=1 Dla n = 1 wezlami kwadratury sa krance przedzialu calkowania, tj. x0 = a, x1 = b. Obliczamy wspolczynniki

A0 = (b a) A1 = (b a)

Z1t 0 Z0 1 t 1

0

1 b a dt = 1 2 b a 0 dt = 0 2

(17)

Kwadratura ta jest zatem rowna wzorowi trapezow (8)

Q(f ) = I (L1 ) =

b a 2

(f (a) + f (b))

(18)

1.3.2 n=2 Dla n = 2 wezlami kwadratury sa x0 = a, x1 = a+2 b , x2 = b. Obliczamy wspolczynniki A0 = A1 =

b a

Z

b a

Z0 2

2 2

A2 = A0

2

0

(t (0 (t (1

Otrzymana kwadratura

Q(f ) = I (L2 ) = jest nazywana wzorem parabol lub

b a 6

1)(t 2) b a dt = 1)(0 2) 6 0)(t 2) 4(b a) dt = 0)(1 2) 6



f (a) + 4f (

wzorem Simpsona.

3

a+b 2

) + f (b)

(19)



(20)

Twierdzenie

Kwadratury Newtona-Cotesa oparte na n + 1 wezlach sa rzedu n + 2 dla n parzystych n + 1 dla n nieparzystych

1.3.3 Zlozone kwadratury Newtona-Cotesa Okazuje sie, ze kwadratury Newtona-Cotesa wyzszych rzedow sa malo przydatne w praktycznych obliczeniach. Na ogol bardziej celowym jest podzielenie przedzialu calkowania [a; b] na N rownych podprzedzialow [xi ; xi+1 ] dlugosci h = bNa punktami xi = a + ih dla i = 0; 1; : : : ; N i stosowanie na kazdym z nich kwadratury Newtona-Cotesa niskiego rzedu. Konstruowane w ten sposob kwadratury, okreslone na calym przedziale [a; b], sa nazywane zlozonymi kwadraturami Newtona-Cotesa.

Jako pierwszy przyklad w kazdym z przedzialow [xi ; xi+1 ] zastosujmy kwadrature Newtona-Cotesa rzedu dwa, tzn. wzor trapezow. Otrzymujemy:

h

TN (f ) = (f (a) + 2

NX1

2

i=1

f (a + ih) + f (b))

(21)

Jako drugi przyklad wyprowadzmy wzor na zlozona kwadrature Simpsona

Z

xi+1

6

xi

Stad

Z a

2

b

h x +x f (x)dx = (f (xi ) + 4f ( i i+1 ) + f (xi+1 )) 2

NX1 NX1 h h f (a + ih) + 4 f (a + (2i 1) ) + f (b)) f (x)dx = (f (a) + 2

6

i=1

i=1

2

(22)

(23)

Problem zbieznosci ciagu kwadratur

Rozpatrzmy ciag kwadratur Qn (n = 0; 1; : : :)

Qn (f ) =

Rb

n X i=0

A(in) f (x(in) )

(24)

przyblizajacych calke Ip (f ) = a p(x)f (x)dx. Interesowac nas bedzie warunki zbieznosci ustalonego ciagu kwadratur do calek funkcji ciaglych na [a; b]. Zakladamy zatem, ze dane sa nieskonczene macierze trojkatne wezlow x(i j ) i wspolczynnikow A(i j )

x(0) 0 (1) x(1) 0 ; x1 (2) (2) x(2) 0 ; x1 ; x2 ::::::::::::

A(0) 0 (1) A0 ; A(1) 1 (2) (2) A0 ; A1 ; A(2) 2 ::::::::::::

de niujace ciag (24). Podstawowym twierdzeniem dajacym odpowiedz na postawione powyzej pytanie jest

Twierdzenie

Ciag kwadratur (24) jest zbiezny dla dowolnych funkcji f ciaglych na odcinku [a; b] lim

n!1

n X i=0

A(in) f (x(in) ) =

wtedy i tyko wtedy, gdy 4

Z

a

b

p(x)f (x)dx

(25)

 

ciag (24) jest zbiezny dla dowolnego wielomianu oraz istnieje taka stala K , ze dla n = 0; 1; : : : zachodzi nierownosc n X i=0

jAin j  K ( )

(26)

Warunki twierdzenia sa spelnione przez zlozone kwadratury Newtona-Cotesa oraz kwadratury Gaussa.

3

Przyspieszanie szybkosci zbieznosci ciagu kwadratur

Zalozmy, ze dla dowolnego n kwadratury Qn spelniaja rownanie

c c 1 Ip (f ) = Qn (f ) + 11 + : : : + mm + O( m+1 ); n n n

(27)

gdzie wykladniki i tworza ciag rosnacy 0 < 1 < 2 < 3 < : : :

(28)

i sa znane, natomiast stale ci zalezne od funkcji f , ale niezalezne od n , sa niewiadome. Idea ekstrapolacji Richardsona polega na wyznaczeniu takiej kombinacji liniowej kwadratur Qn i Qs , ktorej reszta nie zawiera skladnika nc11 . Podstawiajac we wzorze (27) n = ns, przy ustalonym s > 1, dostajemy zaleznosc

c c 1 Ip (f ) = Qsn (f ) + 1 1 1 + : : : + mmm + O( m+1 ) s n s n n Nastepnie mnozac ja stronami przez s1 , odejmujac rownosc (27) i przeksztalcajac otrzymujemy s1 Qsn (f ) Qn (f ) c02 c0 c0 1 Ip (f ) = + 2 + 33 + : : : + mm + O( m+1 )  1 s 1 n n n n

(29)

(30)

gdzie

s1 i 1 c c0i = 1 s 1 i

De niujac

Q1n (f ) =

(31)

s1 Qsn (f ) Qn(f ) s1 1

(32)

okreslamy nowy ciag kwadratur, ktorego reszty sa rzedu n1 2 , a nie n11 jak w przypadku poczatkowego ciagu fQn g. Postepowanie to mazna kontynuowac biorac jako Q2n (f ) odpowiednia kombinacje liniowa Q1n (f ) i Q2sn (f ), ktorej reszta nie zawiera skladnika nc22 , itd. az do Qm n (f ). 0

W przypadku kwadratury trapezow wzor (27) przyjmuje postac

Z

a

b

f (x)dx = Tn(f ) +

c 1 c1 + : : : + 2mm + O( 2m+2 ) n2 n n

(33)

Stosujac do ciagu Tn (f ) ekstrapolacje Richardsona z s = 2 otrzymujemy ciag kwadratur

Tn1(f ) =

4T2n (f ) Tn (f ) 4 1

(34)

Ekstrapolujac k krotnie dostajemy

Tnk =

4k T2kn 1 Tnk 4k 1

5

1

;

Tn0 = Tn

(35)

T10 T20 T40 T80 0 T16 :::

T11 T21 T41 T81 :::

T12 T22 T13 T42 T23 T14 ::: ::: :::

(36)

Kwadratury Tnk nazywane sa kwadraturami Romberga a przedstawiony tu szczegolny przypadek ekstrapolacji Richardsona - metoda Romberga. Latwo sprawdzic, ze kwadratury Tn1 tworzace druga kolumne macierzy (36) sa zlozonymi wzorami Simpsona. Natomiast elementy dalszych kolumn nie sa identyczne z odpowiednimi zlozonymi kwadraturami Newtona-Cotesa.

6