Aproksimacija funkcije i numerička integracija

Lagrangeov interpolacijski polinom

Neka je dana funkcija $f:[a,b]\rightarrow{}R,$ i međusobno različite točke $x_0,x_1,x_2,\ldots,x_n$ u $[a,b].$ Želimo aproksimirati funkciju $f$ polinomom koji u izabranim točkama ima iste vrijednosti kao funkcija $f.$ Ovako postavljen problem ima mnogo rješenja. Međutim, ako zahtijevamo da stupanj polinoma bude najviše $n,$ onda imamo samo jedno rješenje, što upravo tvrdi sljedeći teorem.

Teorem 27   Neka je dana $n+1$ vrijednost

$$f(x_k) = y_k,\hspace{1cm}k=0,1,2,3,\ldots,n$$

funkcije $f:[a,b]\rightarrow{}R.$ Tada postoji jedan i samo jedan polinom $P$ stupnja najviše $n,$ takav da je

$$P(x_k) = y_k,\hspace{1cm}k=0,1,2,3,\ldots,n.$$

Dokaz. Neka su polinomi $L_k, k=0,1,2,\ldots,n$ stupnja najviše $n$ takvi da je

$$L_k(x_j) = \delta_{k\,j},\hspace{1cm}k,j=0,1,2,\ldots,n,$$

gdje je $\delta_{k\,j}$ Kroneckerov simbol 1.1. Tada je

$$P(x) = y_0\,L_0(x) + y_1\,L_1(x) + \cdots + y_n\,L_n(x) = \sum_{i=0}^n y_i\,L_i(x)$$

traženi polinom. Doista,

$$P(x_k) = \sum_{i=0}^n y_i\,L_i(x_k) = \sum_{i=0}^n y_i\,\delta_{i\,k} = y_k.$$

Konstruirajmo sada polinome $L_k.$ Polinom $L_k$ ima nultočke $x_0,x_1,\ldots,x_{k-1},x_{k+1},\ldots,x_n.$ Polinom minimalnog stupnja s tim svojstvom ima oblik

$$L_k(x) = a_k\,(x-x_0)(x-x_1)\cdots (x-x_{k-1})(x-x_{k+1})\cdots (x-x_n).$$

Broj $a_k$ je neodređen. Njega određujemo iz uvjeta $L_k(x_k)=1,$ tj.

$$a_k\,(x_k-x_0)(x_k-x_1)\cdots (x_k-x_{k-1})(x_k-x_{k+1})\cdots (x_k-x_n) = 1.$$

Tako je

$$a_k = \frac{1}{(x_k-x_0)(x_k-x_1)\cdots (x_k-x_{k-1})(x_k-x_{k+1})\cdots (x_k-x_n)},$$

pa je

$$L_k(x) = \frac{(x-x_0)(x-x_1)\cdots ... ...)\cdots (x_k-x_n)} = \prod^n_{\substack{i=0\\ i\neq k}} \frac{x-x_i}{x_k-x_i}.$$

Dakle

$$P(x) = \sum_{k=0}^n y_k\,\prod^n_{\substack{i=0\\ i\neq k}} \frac{x-x_i}{x_k-x_i}.$$

Jedinstvenost ovog polinoma slijedi ovako. Pretpostavimo da je $Q$ također polinom stupnja najviše $n$ takav da je $Q(x_k)=y_k,k=0,1,2,\ldots,n.$ Tada je $f = P-Q$ polinom stupnja najviše $n,$ i ima $n+1$ nultočku. Jedna od posljedica Osnovnog teorema algebre tvrdi da je tada $f=0,$ tj. $f(x)=0,$ za svaki $x\in R.$ Dakle $Q(x)=P(x),$ za svaki $x\in R,$ tj. $P=Q.$ $\heartsuit$

Polinom

$$P(x) = \sum_{k=0}^n f(x_k)\,\prod^n_{\substack{i=0\\ i\neq k}} \frac{x-x_i}{x_k-x_i}$$ (3.13)

se zove Lagrangeov interpolacijski polinom funkcije $f$ za točke $x_0,x_1,x_2,$ $\ldots,x_n.$

Ocjena greške

Radi jednostavnosti pretpostavimo da je

$$x_0<x_1<x_2<\ldots <x_{n-1}<x_n.$$

Teorem 28   Neka funkcija $f$ osim uvjeta iz teorema 27 zadovoljava još uvjet da ima $n+1$-vu neprekidnu derivaciju. Neka je $x\in [a,b]$ proizvoljan, i neka je $S$ najmanji segment koji sadrži točke $x,x_0,x_1,x_2,\ldots,x_n.$ Tada postoji $\xi_x\in S$ takav da je

$$f(x) - P(x) = \frac{f^{(n+1)}(\xi_x)}{(n+1)!}\,(x-x_0)(x-x_1)(x-x_2)\cdots (x-x_n).$$ (3.14)

Dokaz. Ako je $x=x_k$ za neki $k,$ onda nemamo što dokazivati, jer su tada obje strane u (3.14) jednake nuli. Zato pretpostavimo da je $x\in{}[a,b],$ i $x\neq{}x_k$ za svaki $k.$ Stavimo, radi kraćeg zapisa

$$L(x) = (x-x_0)(x-x_1)(x-x_2)\cdots (x-x_n).$$

Formirajmo pomoćnu funkciju

$$F(t) = f(t) - P(t) - c\,L(t),$$ (3.15)

gdje je

$$c = \frac{f(x)-P(x)}{L(x)}.$$

Imamo

$$F(x_k) = f(x_k) - P(x_k) - c\,L(x_k) = 0,\hspace{1cm}k=0,1,2,\ldots,n.$$

Također je

$$F(x) = f(x) - P(x) - c\,L(x) = f(x) - P(x) - \frac{f(x)-P(x)}{L(x)}\,L(x) = 0.$$

Tako funkcija $F$ ima $n+2$ međusobno različite nultočke u $S.$ Prema Rolleovom teoremu njezina derivacija $F'$ ima barem $n+1$ nultočku u $S,$ druga derivacija $F''$ mora imati barem $n$ nultočaka u $S,$ $\ldots,$ $n+1$-va derivacija $F^{(n+1)}$ ima barem jednu nultočku u $S.$ Neka je $\xi_x$ jedna od tih točaka, u kojima se $F^{(n+1)}$ poništava. Derivirajmo (3.15) (po $t$ jer je $t$ varijabla, dok $x$ smatramo fiksnim) $n+1$ puta i uvrstimo $t=\xi_x.$ Budući da je polinom $P$ stupnja najviše $n,$ slijedi $P^{(n+1)}=0.$ Zatim, $L$ je polinom $n+1$-vog stupnja, koeficijent uz $x^{n+1}$ je $1,$ pa je $L^{(n+1)}(t)=(n+1)!.$ Tako imamo

$$0 = F^{(n+1)}(\xi_x) = f^{(n+1)}(\xi_x) - c\,(n+1)!,$$

odnosno

$$c = \frac{f(x)-P(x)}{L(x)} = \frac{f^{(n+1)}(\xi_x)}{(n+1)!}.$$

Odatle

$$f(x)-P(x) = \frac{f^{(n+1)}(\xi_x)}{(n+1)!}\,L(x).$$

$\heartsuit$

Iz ovog teorema proizlazi sljedeća ocjena greške. Neka je

$$M_n=\max_{x\in [a,b]}\left\vert f^{(n+1)}(x)\right\vert.$$

Tada za svaki $x\in [a,b]$ vrijedi

$$\vert f(x)-P(x)\vert \leqslant \frac{M_n}{(n+1)!}\,\vert x-x_0\vert\vert x-x_1\vert\vert x-x_2\vert\cdots \vert x-x_n\vert.$$ (3.16)

Primjer 3.14   Naći Lagrangeov interpolacioni polinom za funkciju $f(x)=e^{-x},$ uzimajući da je

$$e^{-1} = 0.367879,\quad e^{-2} = 0.135335,\quad e^{-3} = 0.0497871,$$

zatim pomoću njega naći približnu vrijednost $f(2.5),$ i ocijeniti grešku.

Rješenje. Po formuli (3.13) imamo

$$P(x) = 0.367879\,\frac{(x-2)(x-3)}{(1-2)(1-3)} + 0.135335\,\frac{(x-3)(x-1)}{(2-3)(2-1)} + 0.0497871\,\frac{(x-1)(x-2)}{(3-1)(3-2)}$$

$$= 0.747419 - 0.453038\,x + 0.073498\,{x^2}.$$

Odatle

$$f(2.5) = e^{-2.5} = 0.0741865.$$

Ocijenimo sada grešku.

$$M_2=\max_{x\in [1,3]}\left\vert f^{(3)}(x)\right\vert = \max_{x\in [1,3]} e^{-x} = e^{-1} = 0.367879.$$

Dakle, prema formuli (3.16), greška koju pri tom činimo nije veća od

$$\frac{0.36788}{6}\,\vert 2.5-1\vert\vert 2.5-2\vert\vert 2.5-3\vert = 0.0229925.$$

Metoda najmanjih kvadrata

Na sljedećoj slici je funkcija, zadana svojim vrijednostima u nekoliko točaka, interpolirana polinomom prvog stupnja (pravac) i polinomom drugog stupnja (crtkana linija) u smislu metode najmanjih kvadrata.

Lagrangeov polinom veoma dobro aproksimira funkciju lokalno, u izabranim točkama, dok izvan tih točaka aproksimacija može biti vrlo loša. Sada ćemo upoznati metodu najmanjih kvadrata, metodu pomoću koje možemo zadanu funkciju aproksimirati drugom funkcijom određenog tipa globalno, tako da u izvjesnom smislu njihova međusobna udaljenost bude što manja, bez obzira na to što se funkcije možda neće poklapati niti u jednoj točki.

Pretpostavimo najprije da su vrijednosti funkcije poznate samo u nekim točkama. Neka su $x_1, x_2, \ldots, x_n$ dane točke i neka su $f(x_1),f(x_2),\ldots{},f(x_n)$ pripadne vrijednosti funkcije $f.$ Želimo naći onu funkciju $g$ određenog tipa s neodređenim parametrima $\alpha,\beta,\ldots{},$ koja najbolje aproksimira funkciju $f.$ To možemo učiniti na sljedeći način.

Izračunamo sumu kvadrata razlika funkcija $f$ i $g$

$$S = \sum_{k=1}^n \left[f(x_k) - g(x_k)\right]^2,$$

i zatim tražimo neodređene parametre $\alpha,\beta,\ldots{}$ iz uvjeta da $S$ kao funkcija tih parametara ima minimalnu vrijednost.

Ako je funkcija $f$ zadana u svim točkama nekog segmenta $[a,b],$ onda se funkcija $S$ definira pomoću integrala

$$S = \int_{a}^b \left[f(x) - g(x)\right]^2\,dx.$$

Klase funkcija iz kojih biramo funkciju $g$ su obično polinomi prvog stupnja $g(x)=\alpha{}x+\beta{},$ polinomi drugog stupnja $g(x)=\alpha{}x^2 +\beta{}x+\gamma{},$ eksponencijalne funkcije $g(x)=\alpha{}\,e^{\beta{}\,x}+\gamma{},$ itd.

Na pr. ako želimo naći polinom prvog stupnja najbliži funkciji $f$ čije su nam vrijednosti poznate u točkama $x_1,x_2,\ldots{},x_n,$ onda je

$$S(\alpha{},\beta{}) = \sum_{k=1}^n \left[f(x_k) - \alpha{}\,x_k - \beta{}\right]^2.$$ (3.17)

Nužan uvjet za ekstrem funkcije $S$ je

$$\frac{\partial{}S}{\partial{}\alpha{}} = 0,\hspace{1cm} \frac{\partial{}S}{\partial{}\beta{}} = 0,$$

dakle

$$\sum_{k=1}^n \left[f(x_k) - \alpha{}\,x_k - \beta{}\right]\,x_k = 0,$$

$$\sum_{k=1}^n \left[f(x_k) - \alpha{}\,x_k - \beta{}\right] = 0.$$

Kad sredimo ovaj sustav jednadžbi, dobijemo

$$\alpha{}\,\sum_{k=1}^n x_k^2 + \beta{}\,\sum_{k=1}^n x_k = \sum_{k=1}^n f(x_k)\,x_k,$$

$$\alpha{}\,\sum_{k=1}^n x_k + n\,\beta{} = \sum_{k=1}^n f(x_k).$$

Iz njega izračunamo $\alpha$ i $\beta.$ Iz prirode problema jasno je da funkcija $S$ ima samo minimum, pa tako određeni parametri doista daju najbolju aproksimaciju funkcioje $f.$

Kad se funkcija $g$ traži u obliku eksponencijalne, logaritamske i slično, onda se na ovaj način u pravilu dobije sustav nelinearnih jednadžbi. Takve sustave je teže rješavati nego linearne. Obično se nekim operacijama nad funkcijama (na pr. logaritmiranjem) pojednostavni klasa funkcija u kojoj tražimo aproksimaciju. No, parametri koje tako dobijemo nisu najbolji mogući u smislu metode najmanjih kvadrata.

Primjer 3.15   Neka su funkcija i točke kao u primjeru 3.14. Interpolirati funkciju polinomom drugog stupnja, ali metodom najmanjih kvadrata, i pomoću njega približno naći $f(2.5).$

Rješenje. Neka je

$$p(x) = a + b\,x + c\,x^2$$

polinom koji u smislu metode najmanjih kvadrata interpolira podatke

$x$	$1$	$2$	$3$
$f(x)$	$0.367879$	$0.135335$	$0.0497871$

Sada formulu (3.17) treba napisati u obliku

$$S(a,b,c) = \sum_{k=1}^3 \left[f(x_k) - c\,x^2 - b\,x_k - a\right]^2.$$

Ova funkcija ovisi o tri varijable, pa je nužan uvjet za ekstrem

$$\frac{\partial S(a,b,c)}{\partial a} = -2\,\sum_{k=1}^3 \left[f(x_k) - c\,x^2 - b\,x_k - a\right] = 0$$

$$\frac{\partial S(a,b,c)}{\partial b} = -2\,\sum_{k=1}^3 \left[f(x_k) - c\,x^2 - b\,x_k - a\right]\,x_k = 0$$

$$\frac{\partial S(a,b,c)}{\partial c} = -2\,\sum_{k=1}^3 \left[f(x_k) - c\,x^2 - b\,x_k - a\right]\,x_k^2 = 0.$$

Kad uvrstimo poznate veličine i sredimo, dobivamo sustav jednadžbi

$$6\,a + 12\,b + 28\,c = 1.106$$

$$12\,a + 28\,b + 72\,c = 1.57582$$

$$28\,a + 72\,b + 196\,c = 2.71461$$

Rješenje je $a=0.74742,b=-0.453038,c=0.073498,$ pa je traženi polinom

$$p(x) = 0.74742 - 0.453038\,x + 0.073498\,{x^2}.$$

Tražena vrijednost funkcije je

$$f(2.5) = e^{-2.5} = 0.0741865.$$

Numerička integracija

Zadatak je izračunati integral

$$\int_a^b f(x)\,dx.$$

Umjesto da integriramo podintegralnu funkciju, što često nije moguće, ili zahtijeva puno posla, integriramo polinom, koji interpolira funkciju u odgovarajućim točkama. Numerički to radimo na sljedeći način. Segment $[a,b]$ podijelimo na podsegmente točkama

$$a=x_0<x_1<x_2<\ldots<x_{n-1}<x_n=b.$$

Radi jednostavnosti i određenosti postupka podjela se uzme ekvidistantnom. U tim točkama interpoliramo funkciju Lagrangeovim polinomom, i zatim integriramo polinom. Tako dobiven broj predstavlja približnu vrijednost zadanog integrala. Dakle

$$\int_a^b f(x)\,dx \approx{} \int_a^b P(x)\,dx.$$

Ako uvrstimo Lagrangeov polinom (3.13), imamo

$$\int_a^b f(x)\,dx \approx{} \int_a^b... ...f(x_k)\,\int_a^b \prod^n_{\substack{i=0\\ i\neq k}} \frac{x-x_i}{x_k-x_i}\,dx.$$ (3.18)

Da bismo razmatranje učinili neovisnim o segmentu $[a,b],$ svedimo ga supstitucijom na fiksni segment $[-1,1].$ To možemo učiniti afinom funkcijom (polinomom prvog stupnja) čiji je graf pravac kroz točke $(a,-1)$ i $(b,1).$

Jednadžba tog pravca je

$$t = -1 + \frac{2}{b-a}\,(x-a),$$

odnosno

$$x = \frac{b-a}{2}\,t + \frac{b+a}{2}.$$

Supstitucija čuva ekvidistantnost, pa je

$$-1=t_0<t_1<t_2<\ldots,t_{n-1}<t_n=1$$

ekvidistantna podjela segmenta $[-1,1]$ na $n$ podsegmenata. Tom supstitucijom formula (3.18) prelazi u

$$\int_a^b f(x)\,dx \approx{} \frac{b-a}{2}\,\sum_{k=0}^n w_k\, f\left(\frac{b-a}{2}\,t_k + \frac{b+a}{2}\right),$$

gdje je

$$w_k = \int_{-1}^1 \prod^n_{\substack{i=0\\ i\neq k}} \frac{t-t_i}{t_k-t_i}\,dt,\hspace{1cm}k=0,1,2,\ldots,n.$$ (3.19)

Vidimo da ponderi $w_k$ ne ovise o segmentu, niti o funkciji, već samo o broju $n.$

Iz formule (3.14) slijedi da je greška koju pri tom činimo

$$R_n = \int_a^b f(x)\,dx - \int_a^b P(x)\,dx = \frac{f^{(n+1)}(\xi_x)}{(n+1)!}\,\int_a^b L(x)\,dx,$$

gdje je

$$L(x) = (x-x_0)(x-x_1)(x-x_2)\cdots (x-x_n).$$

Točka $\xi_x$ nam nije poznata, pa za ocjenu greške moramo uzeti maksimum ove derivacije na segmentu $[a,b]$

$$\vert R_n\vert \leqslant{} \frac{M_n}{(n+1)!}\,\left\vert\int_a^b... ...t\vert \leqslant{} \frac{M_n}{(n+1)!}\,\int_a^b \left\vert L(x)\right\vert\,dx,$$

gdje je

$$M_n=\max_{x\in [a,b]}\left\vert f^{(n+1)}(x)\right\vert.$$

Specijalno za $n=1$ imamo $x_0=a<x_1=b.$ Zatim,

$$L(x) = (x-a)(x-b),$$

pa je

$$\vert R_1\vert \leqslant{} \frac{M_1}{(n+1)!}\int_a^b (x-a)(b-x)\,dx = \frac{M_1\,\left(b-a \right)^3}{12},$$

gdje je

$$M_1 = \max_{x\in [a,b]} \left\vert f''(x)\right\vert.$$

Za $n=2,$ imamo tri točke podjele $x_0=a<x_1<x_2=b.$ Pretpostavimo da je podjela ekvidistantna, tj. da je $x_1=\frac{a+b}{2}.$ Tada je

$$L(x) = \frac{1}{2}(x-a)(2\,x-(a+b))(x-b).$$

Kad bismo ponovili postupak kao gore, dobili bismo ocjenu reda veličine $(b-a)^4$ No, ta se ocjena može poboljšati. Zbog činjenice da je funkcija $L(x)$ simetrična u odnosu na točku $\frac{a+b}{2}$ njezin integral po segmentu $[a,b]$ iščezava. U tom slučaju se funkcija može interpolirati s polinomom 4. stupnja, tako da se točka $\frac{a+b}{2}$ uzme kao dvostruka. O tome kako se to radi ovdje nećemo govoriti. Primijetimo samo da zbog integriranja polinoma 4. stupnja ocjena postaje reda veličine $(b-a)^5.$ Može se pokazati de je ona

$$\vert R_2\vert \leqslant{} \frac{M}{90}\,\left(\frac{b-a}{2}\right)^5 = \frac{M\,(b-a)^5}{2880},$$

gdje je

$$M = \max_{x\in [a,b]} \left\vert f^{(iv)}(x)\right\vert.$$

Trapezna formula

U slučaju $n=1,$ iz formule (3.19) imamo $t_0=-1,t_1=1,$ pa su ponderi

$$w_0 = \int_{-1}^1 \frac{t-1}{-2}\,dt = 1,\hspace{1cm}w_1 = \int_{-1}^1 \frac{t+1}{2}\,dt = 1.$$

Prema tome formula glasi

$$\int_a^b f(x)\,dx \approx{} \frac{b-a}{2}\,\left[f(a)+f(b)\right],$$ (3.20)

uz ocjenu greške

$$\vert R_1\vert \leqslant{} \max_{x\in [a,b]} \left\vert f''(x)\right\vert\frac{(b-a)^3}{12}.$$

Općenito je greška ove aproksimacije velika. Zato dijelimo segment $[a,b]$ na podsegmente, i na svakom od njih posebno koristimo ovu formulu. Radi jednostavnosti uzimamo ekvidistantnu podjelu. Neka je dakle

$$\frac{b-a}{m}=h,$$   i $$\hspace{1cm}a=x_0<x_1<x_2<\ldots{}<x_{m-1}<x_m=b,$$

gdje je $x_i=a+ih.$ Ako na svaki od podsegmenata, čija je duljina $h,$ primijenimo formulu (3.20), dobivamo

$$\int_a^b f(x)\,dx \approx{} \frac{h}... ...ft[f(x_1)+f(x_2)\right] + \cdots{} + \frac{h}{2}\,\left[f(x_{m-1})+f(b)\right],$$

odnosno

$$\int_a^b f(x)\,dx \approx{} \frac{h}... ...eft\{f(a)+2\,\left[f(x_1) + f(x_2) + \cdots{} + f(x_{m-1})\right]+f(b)\right\}.$$ (3.21)

Formula (3.21) se zove trapezna formula.

Ocjena greške se dobije zbrajanjem grešaka na pojedinim podsegmentima. Tako imamo

$$\vert R_T\vert \leqslant{} \sum_{j=1}^m \frac{M_j\,h^3}{12},$$

gdje je

$$M_j = \max_{x\in [x_{j-1},x_j]} \left\vert f''(x)\right\vert,\hspace{1cm}j=1,2,\ldots{},m.$$

Za ovu ocjenu trebalo bi mnogo puta ($m$ puta) ocjenjivati drugu derivaciju i to na raznim podsegmentima, što nije jednostavno. Zato radije ocjenu napravimo jednom za cijeli segment $[a,b],$ i uzmemo u obzir da je za svaki $j$

$$M_j \leqslant{} M = \max_{x\in [a,b]} \left\vert f''(x)\right\vert.$$

Tako imamo ocjenu

$$\vert R_T\vert \leqslant{} \sum_{j=1}^m \frac{M\,h^3}{12} = \frac{M\,h^3}{12}\,m = \frac{M\,(b-a)^3}{12\,m^2}.$$

Simpsonova formula

U slučaju $n=2,$ iz formule (3.19) imamo $t_0=-1,t_1=0,t_2=1$ pa su ponderi

$$w_0 = \int_{-1}^1 \frac{t(t-1)}{2}\,dt = \frac{1}{3},\hspace{1cm}w_1 = -\int_{-1}^1 (t+1)(t-1)\,dt = \frac{4}{3},$$

$$w_2 = \int_{-1}^1 \frac{t(t+1)}{2} = \frac{1}{3}.$$

Prema tome formula glasi

$$\int_a^b f(x)\,dx \approx{} \frac{b-a}{6}\,\left[f(a) + 4\, f\left(\frac{a+b}{2}\right) + f(b)\right],$$

uz ocjenu greške

$$\vert R_2\vert \leqslant{} \frac{M}{90}\,\left(\frac{b-a}{2}\right)^5,$$

gdje je

$$M = \max_{x\in [a,b]} \left\vert f^{(iv)}(x)\right\vert.$$

Kao i kod trapezne formule, točniji rezultat ćemo dobiti, ako segment $[a,b]$ podijelimo na podsegmente, i na svakom od njih posebno koristimo ovu formulu. Radi jednostavnosti uzimamo ekvidistantnu podjelu. Budući da na svakom podsegmentu uzimamo još srednju točku, podsegmenata imamo zapravo $2m.$ Neka je dakle

$$\frac{b-a}{2\,m}=h,$$   i $$\hspace{1cm} a=x_0<x_1<x_2<\ldots{}<x_{2m-1}<x_{2m}=b.$$

Sada primjenjujemo približnu formulu na parovima susjednih podsegmenata i to tako da srednja točka bude točka s neparnim indeksom. Tako dobivamo formulu

$$\int_a^b f(x)\,dx \approx{} \frac{b-a}{6\,m}\,\left[f(a) + 4\, f(x_1) + f(x_2)\right]+$$

$$+ \frac{b-a}{6\,m}\,\left[f(x_2) + 4\, f(x_3) + f(x_4)\right] + \cdots{}+$$

$$+ \frac{b-a}{6\,m}\,\left[f(x_{2m-2}) + 4\, f(x_{2m-1}) + f(x_{b})\right],$$

odnosno

$$\int_a^b f(x)\,dx\approx\frac{h}{3}\{f(a)+2\,[f(x_2)+f(x_4)+\cdots +f(x_{2m-2})]$$

$$+4\,[f(x_1)+f(x_3)+\cdots+f(x_{2m-1})]+f(b)\},$$

koja se zove Simpsonova formula.

Diskusija kao kod trapezne formule nas dovodi do ocjene greške

$$\vert R_S\vert \leqslant{} \max_{x\in [a,b]} \left\vert f^{(iv)}(... ...x_{x\in [a,b]} \left\vert f^{(iv)}(x)\right\vert \,\frac{(b-a)^5}{180\,(2m)^4}.$$

Primjer 3.16   Pomoću numeričke integracije izračunati približno broj $\pi$ na šest decimala točno.

Rješenje. Najprije sjetimo se da je ${\rm Arctg}\,1=\frac{\pi}{4},$ i ${\rm Arctg}\,0=0.$ Odatle

$$\frac{\pi}{4} = {\rm Arctg}\,1 - {\r... ...= \int_0^1 \frac{d\,{\rm Arctg}\,x}{d\,x}\,dx = \int_0^1 \frac{1}{1 + x^2}\,dx.$$

Dakle

$$\pi = 4\,\int_0^1 \frac{1}{1 + x^2}\,dx,$$

pa treba izračunati ovaj integral na šest decimala točno. Račun ćemo provesti na dva načina, pomoću trapezne i pomoću Simpsonove formule.

1. Podintegralna funkcija je

$$f(x) = \frac{4}{1 + x^2},$$

njezina druga derivacija je

$$f''(x) = {\frac{8\,\left( -1 + 3\,{x^2} \right) } {{{\left( 1 + {x^2} \right) }^3}}},$$

a treća

$$f'''(x) = {\frac{-96\,x\,\left( -1 + {x^2} \right) }{{{\left( 1 + {x^2} \right) }^4}}}.$$

Kako treća derivacija nema nula u intervalu $\langle 0,1\rangle,$ druga derivacija je na tom intervalu monotona. Računanjem vrijednosti na rubovima lako se vidi da ona raste i to od $-8$ do $2.$ To znači da je $M=8.$ Dakle imamo ocjenu greške za trapeznu formulu

$$\vert R_T\vert \leqslant{} \frac{8}{12\,m^2}.$$

Budući da tražimo točnost prvih šest decimala, mora biti

$$\frac{8}{12\,m^2} < 0.5\times 10^{-6},$$

tj.

$$m > \sqrt{\frac{8\cdot 10^{6}}{12\cdot 0.5}} = 1154.7$$

Prema tome treba segment $[0,1]$ podijeliti na $1155$ podsegmenata.

2. Pomoću Simpsonove formule račun ide ovako. Četvrta derivacija je

$$f^{iv}(x) = {\frac{96\,\left( 1 - 10\,{x^2} + 5\,{x^4} \right) }{{{\left( 1 + {x^2} \right) }^5}}},$$

a peta

$$f^{v} = {\frac{-960\,x\,\left( 3 - 10\,{x^2} + 3\,{x^4} \right) }{{{\left( 1 + {x^2} \right) }^6}}}.$$

Peta derivacija se poništava u $x=0.57735,$ ali je $f^{vi}$ u toj točki pozitivna, pa $f^{iv}$ ima u $x=0.57735$ minimum. Tako četvrta derivacija pada na intervalu od 0 do $0.57735,$ i zatim raste na intervalu od $0.57735$ do $1.$ Imamo

$$f^{iv}(0) = 96,\quad f^{iv}(0.57735) = -40.5,\quad f^{iv}(1) = -12.$$

Dakle $M=96,$ i prema tome treba biti

$$\frac{96}{2880\,m^4} < 0.5\times 10^{-6},$$

odakle $m>16.0686.$ Budući da je broj segmenata potreban za Simpsonovu formulu paran, minimalni $m$ koji treba uzeti je $18.$ Kad se provede potreban račun dobije se

$$\pi \approx 3.141592652.$$

Gaussova kvadratura

Trapezna i Simpsonova formula su koristile ekvidistantnu podjelu. Točke podjele su bile jednoliko razmještene na segmentu. Te formule u pravilu računaju točno integrale polinoma najviše onog stupnja kojeg je bio interpolacioni polinom. Tako trapezna formula računa točno integrale polinoma do uključivo prvog stupnja.

Postavlja se pitanje da li postoje formule koje računaju točno integrale polinoma stupnja višeg nego što je interpolacioni polinom. Takve formule postoje i zovu se Gaussove kvadraturne formule. U tim formulama točke podjele nisu više ekvidistantno raspoređene. Točke podjele ćemo u daljnjem zvati čvorovima.

Kvadraturna formula općenito je oblika

$$\int_{a}^{b}\,f(x)\,dx\approx \frac{b-a}{2}\,\sum_{j=1}^n c_j\,f(x_j).$$ (3.22)

Zbog svojstva linearnosti integrala, kvadraturna formula koja računa točno potencije do uključivo stupnja $m,$ računat će točno i njihove linearne kombinacije, tj. polinome do uključivo $m$-tog stupnja. Dakle, pretpostavimo da greška prilikom računanja potencija do uključivo stupnja $m$ iščezava

$$R_n(x^k)=\int_{a}^{b}\,x^k\,dx- \frac{b-a}{2}\, \sum_{j=1}^n c_j\,x^k_j=0,\hspace{1cm}k=0,1,2,\ldots,m.$$

Tako imamo $m+1$ jednadžbu za $2n$ nepoznanica $c_1,c_2,\ldots,c_n,x_1,x_2,\ldots,x_n.$ Općenito sustav treba imati onoliko jednadžbi koliko ima nepoznanica, da bismo imali jedinstveno rješenje. Ne ulazeći dublje u tu problematiku, primijetimo da u ovom slučaju to znači da treba biti $m\leq 2n-1.$ Dakle, u principu, s $n$ kvadraturnih točaka $x_1, x_2, \ldots, x_n$ se mogu rješavati točno integrali polinoma stupnja najviše $2n-1.$

Umjesto rješavanja gornjeg sustava, za određivanje čvorova $x_1, x_2, \ldots, x_n$ i težina $c_1,c_2,\ldots,c_n$ koristimo ideju koja potječe od C. F. Gaussa, a osniva se na sljedećem razmatranju.

Neprekidne funkcije na $[a,b]$ čine vektorski prostor obzirom na uobičajene operacije zbrajanja funkcija i množenja funkcije brojem. Svakom paru $f,g$ takvih funkcija možemo pridružiti broj

$$f\cdot g=\int_a^b f(x)\,g(x)\,dx.$$

Funkcija, koja uređenom paru $(f,g)$ pridružuje $f\cdot{}g,$ ima svojstva

1.

$(f+g)\cdot h = \int_a^b (f(x)+g(x))\,h(x)\,dx = \int_a^b f(x)\,h(x)\,dx + \int_a^b g(x)\,h(x)\,dx = f\cdot h + g\cdot h,$

2.

$(\lambda\,f)\cdot{}g = \int_a^b \lambda{}f(x)\,g(x)\,dx = \lambda{}\int_a^b f(x)\,g(x)\,dx = \lambda{}f\cdot g,$

3.

$f\cdot g=\int_a^b f(x)\,g(x)\,dx = \int_a^b g(x)\,f(x)\,dx = g\cdot f,$

4.

$f\cdot f=\int_a^b [f(x)]^2\,dx \geqslant{} 0,$

5.

$f\cdot f=\int_a^b [f(x)]^2\,dx = 0\hspace{1cm}\Leftrightarrow \hspace{1cm}f=0.$

pa se zato zove skalarni produkt. To nam omogućava da govorimo o međusobno ortogonalnim funkcijama, njihovoj duljini itd.

Vratimo se sada na naš problem. Da diskusija ne bi ovisila o području integracije, prijeđimo na segment $[-1,1]$ supstitucijom

$$x = \frac{b-a}{2}\,t + \frac{b+a}{2}.$$

Tada je

$$\int_{a}^{b}\,f(x)\,dx = \frac{b-a}{2}\,\int_{-1}^{1}\,f\left(\fr... ...}\,t + \frac{b+a}{2}\right)\,dt = \frac{b-a}{2}\,\int_{-1}^{1}\,\varphi(t)\,dt,$$

gdje je

$$\varphi(t) = f\left(\frac{b-a}{2}\,t + \frac{b+a}{2}\right) = f(x).$$

Iz formule (3.22) slijedi

$$\int_{-1}^{1}\,\varphi(t)\,dt \appro... ...eft(\frac{b-a}{2}\,t_j + \frac{b+a}{2}\right) = \sum_{j=1}^n c_j\,\varphi(t_j).$$

Neka su $t_1,t_2,\ldots,t_n$ tražene kvadraturne točke. Po pretpostavci, kvadraturna formula je točna za polinome do uključivo $2n-1$-og stupnja. Prema tome

$$\int_{-1}^1 (t-t_1)(t-t_2)\cdots(t-t_n)t^k\,dt- \sum_{j=1}^n c_j(t_j-t_1)(t_j-t_2)\cdots(t_j-t_n)t_j^k=0,$$

za $k=0,1,2,\ldots,n-1.$ S druge strane suma iščezava, jer se za svaki $j$ jedan faktor poništava. To znači, da je polinom

$$P_n(t)=(t-t_1)(t-t_2)\cdots(t-t_n)$$

okomit na polinome $1,t,t^2,\ldots,t^{n-1},$ a prema tome i na njihove linearne kombinacije, tj. na sve polinome stupnja najviše $n-1.$ Takvi polinomi se zovu Legendreovi polinomi. Točke Gaussove kvadrature su upravo korijeni (nule) tih polinoma. Nađimo Legendreove polinome. U tu svrhu ortogonalizirajmo polinome $1,t,t^2,\ldots$ Gram-Schmidtovim postupkom

$$P_0(t) = 1,$$

$$P_1(t) = t - \frac{t\cdot{}P_0}{P_0\cdot{}P_0}\,P_0 = t - \frac{\int_{-1}^{1}\, t\,dt}{\int_{-1}^{1}\,dt}= t,$$

$$P_2(t) = t^2 - \frac{t^2\cdot{}P_0}{P_0\cdot{}P_0}\,P_0 - \frac{t... ...\frac{\int_{-1}^{1}\,t^3\,dt}{\int_{-1}^{1}\,t^2\,dt}\,t = {t^2}-{\frac{1}{3}},$$

$$P_3(t) = t^3 - \frac{t^3\cdot{}P_0}{P_0\cdot{}P_0}\,P_0 - \frac{t... ...\frac{t^3\cdot{}P_2}{P_2\cdot{}P_2}\,P_2 = \ldots{} = {t^3} - {\frac{3\,t}{5}},$$

$$P_4(t) = \ldots{} = {t^4} - {\frac{6\,{t^2}}{7}} + {\frac{3}{35}},$$

$$\ldots$$

Tako na segmentu $[-1,1]$ imamo točke kvadrature

$n$	$t_1$	$t_2$	$t_3$	$t_4$
$1$	$$0$$
$2$	$-0.57735$	$0.57735$
$3$	$-0.774597$	$$0$$	$0.774597$
$4$	$-0.861136$	$-0.339981$	$0.339981$	$0.861136$

Težine $c_j$ jednostavno možemo izračunati na sljedeći način. Gaussova kvadratura računa točno polinome do uključivo $2n-1$-og stupnja. Neka su $L_k^{(n)}, k=1,2,3,\ldots,n$ polinomi sa svojstvom

$$L_k^{(n)}(t_j) = \delta_{k\,j},\hspace{1cm}k,j=1,2,\ldots,n.$$

To su polinomi koji su se pojavili prilikom Lagrangeove interpolacije. Njihov stupanj je $n-1,$ pa Gaussova kvadratura računa njihove integrale točno. To znači

$$\int_{-1}^{1}\,L_k^{(n)}(t)\,dt = \sum_{j=1}^n c_j\,L_k^{(n)}(t_j) = c_k,\hspace{1cm} k=1,2,\ldots,n,$$

tj.

$$c_k = \int_{-1}^{1}\,L_k^{(n)}(t)\,dt,\hspace{1cm} k=1,2,\ldots,n.$$

$n$	$c_1$	$c_2$	$c_3$	$c_4$
$1$	$2.000000$
$2$	$1.000000$	$1.000000$
$3$	$0.555556$	$0.888889$	$0.555556$
$4$	$0.347855$	$0.652145$	$0.652145$	$0.347855$

Nultočke Legendreovihovih polinoma uglavnom su iracionalni brojevi, pa uzimajući vrijednosti iz gornje tablice, činimo grešku zamjenjujući iracionalan broj decimalnim. Tako ipak ne dobivamo apsolutnu točnost. No, barem smo uklonili grešku koja se odnosi na zamjenu integrala kvadraturnom formulom.