Rješavanje parcijalnih diferencijalnih jednadžbi

Metoda konačnih diferencija

Rubne i rubno-početne probleme, u kojima se pojavljuju parcijalne diferencijalne jednadžbe, možemo također rješavati metodom konačnih diferencija, ili kako se u ovom slučaju često kaže metodom mreže. Osnovna ideja je zamijeniti derivacije (u ovom slučaju parcijalne) s konačnim diferencijama, i time svesti parcijalnu diferencijalnu jednadžbu na sustav algebarskih jednadžbi. Pokažimo kako se to u pojedinim slučajevima radi.

Jednadžba ravnoteže

Neka je dan rubni problem

$$\Delta\,u(x,y) = f(x,y),\quad \Omega{}=[a,b]\times{}[c,d]$$

$$u(x,y) = g(x,y),\quad \partial{}\Omega{},$$

gdje je $f$ neprekidna funkcija na $\Omega{},$ a $g$ neprekidna na rubu $\partial{\Omega{}}.$

Područje $\Omega$ prekrijemo mrežom koju dobijemo tako da segmente na osi $x$ i osi $y$ podijelimo ekvidistantno na podsegmenate. Tako imamo

$$h = \frac{b-a}{n},\hspace{1cm}k = \frac{d-c}{m}.$$

Broj $h$ se zove korak mreže po $x$ osi, a $k$ korak mreže po $y$ osi. Tako imamo točke podjele

$$x_i = a + i\,h,\quad i=0,1,2,\ldots{},n\hspace{1cm}y_j = c + j\,k,\quad j=0,1,2,\ldots{},m$$

na osima. Točku s koordinatama

$$(x_i,y_j) = (a + i\,h,c + j\,k)$$

zovemo $ij$-tim čvorom mreže. Čvor zovemo unutrašnjim, ako je $1\leqslant{}i\leqslant{}n-1$ i $1\leqslant{}j\leqslant{}m-1.$ U protivnom kažemo da je čvor rubni. Stavimo

$$u(x_i,y_j) = u_{i\,j},\hspace{1cm}f(x_i,y_j) = f_{i\,j}.$$

Kao i u slučaju običnih derivacija, aproksimacije parcijalnih derivacija možemo dobiti pomoću Taylorove formule za funkcije od dvije varijable. Tako imamo

$$\frac{\partial{}^2u(x_i,y_j)}{\parti... ..._j)}{\partial{}y^2} \approx{} \frac{u_{i\,j-1} -2\,u_{i\,j} +u_{i\,j+1}}{k^2}.$$

Neka je $ij$ unutrašnji čvor. Tada diferencijalnu jednadžbu u tom čvoru možemo zamijeniti algebarskom jednadžbom

$$\frac{u_{i-1\,j} -2\,u_{i\,j} +u_{i+1\,j}}{h^2} + \frac{u_{i\,j-1} -2\,u_{i\,j} +u_{i\,j+1}}{k^2} = f_{i\,j}.$$

Ako je čvor $ij$ rubni, onda imamo

$$u_{i\,j} = g_{i\,j}.$$

Time smo dobili onoliko linearnih algebarskih jednadžbi koliko imamo nepoznanica $u_{i\,j}.$ Rješavati treba samo sustav jednadžbi unutrašnjih čvorova. Pretpostavimo, radi jednostavnosti, da su $a,b,c,d,n,m$ takvi da je $h=k.$ Tada, nakon množenja s $h^2,$ jednadžba $ij$-tog čvora postaje

$$u_{i\,j-1}+u_{i-1\,j} -4\,u_{i\,j} +u_{i+1\,j} + u_{i\,j+1} = h^2\,f_{i\,j}.$$

Vidimo da u svakom retku matrice ima najviše pet elemenata različitih od nule. Na sljedećoj slici su istaknuti oni čvorovi, koji imaju svojstvo da vrijednosti funkcije $u$ u njima nalaze u jednadžbi.

Jednadžbe trebamo na neki način poredati. To možemo učiniti tako da čvorove poredamo. Njih možemo poredati na različite načine. Jedan od njih je poredak kao na slici

Dakle

$$1\,1,\quad 2\,1,\quad 1\,2,\quad 3\,1,\quad 2\,2,\quad 1\,3,\quad4\,1,\quad3\,2,\quad2\,3,\quad1\,4,\quad\ldots{}\ .$$

Stavimo

$$u_{\text{\footnotesize\ding{192}}} =... ... = u_{1\,2},\quad u_{\text{\footnotesize\ding{195}}} = u_{3\,1},\quad \ldots\ .$$

Na taj način dobivamo sljedeći sustav linearnih algebarskih jednadžbi

$$-4\,u_{\text{\footnotesize\ding{192}}} + u_{\text{\footnotesize\ding{193}}} + u_{\text{\footnotesize\ding{194}}} = h^2\,f_{1\,1} - g_{1\,0} - g_{0\,1}$$

$$u_{\text{\footnotesize\ding{192}}} -... ...93}}} + u_{\text{\footnotesize\ding{195}}} + u_{\text{\footnotesize\ding{196}}} = h^2\,f_{2\,1} - g_{2\,0}$$

$$u_{\text{\footnotesize\ding{192}}} -... ...94}}} + u_{\text{\footnotesize\ding{196}}} + u_{\text{\footnotesize\ding{197}}} = h^2\,f_{1\,2} - g_{0\,2}$$

$$u_{\text{\footnotesize\ding{193}}} -... ...95}}} + u_{\text{\footnotesize\ding{198}}} + u_{\text{\footnotesize\ding{199}}} = h^2\,f_{3\,1} - g_{3\,0}$$

$$u_{\text{\footnotesize\ding{193}}} +... ...96}}} + u_{\text{\footnotesize\ding{199}}} + u_{\text{\footnotesize\ding{200}}} = h^2\,f_{2\,2}$$

$$\ldots$$

Prilikom izbora poretka treba paziti na to da matrica dobivenog sustava jednadžbi bude što uža.

Jednadžba provođenja

Jednodimenzionalno provođenje topline opisano je sljedećim rubno-početnim problemom

$$\frac{\textstyle{\partial u(x,t)}}{\textstyle{\partial t}} = c^2\,\frac{\textstyle{\partial^2 u(x,t)}}{\textstyle{\partial x^2}},\quad x\in [0,\ell],\;t\in [0,\infty{}\rangle{},$$

$$u(0,t)=0,\quad u(\ell,t)=0,\quad t\in [0,\infty{}\rangle{},$$

$$u(x,0) = g(x),\quad x\in [0,\ell].$$

Neka je

$$h = \frac{\ell}{n}$$

korak mreže po osi $x,$ a $\tau$ korak po osi $t.$ Tako imamo točke podjele na $x$ osi, odnosno $t$ osi

$$x_i = i\,h,\quad i=0,1,2,\ldots{},n\hspace{1cm}t_j = j\,\tau,\quad j=0,1,2,\ldots{}.$$

Interesira nas što se događa u čvorovima, tj. točkama

$$(x_i,t_j) = (i\,h,j\,\tau),\qquad i=0,1,2,\ldots{},n,\;j=0,1,2,\ldots{}\ .$$

Stavimo

$$u(x_i,t_j) = u_{i\,j}.$$

Iz rubnih i početnih uvjeta slijedi

$$u_{0\,j} = u_{\ell\,j} = 0,\quad j=0,1,2,\ldots{},\qquad u_{i\,0} = g(x_i) = g_i,\quad i=0,1,2,\ldots{},n.$$

Dakle ostaje odrediti $u_{i\,j}$ u unutrašnjim čvorovima. U tu svrhu derivacije aproksimiramo koristeći neke od mogućih formula, na pr.

$$\frac{\textstyle{\partial u(x_i,t_j)... ..._j)}{\partial{}x^2} \approx{} \frac{u_{i-1\,j} -2\,u_{i\,j} +u_{i+1\,j}}{h^2}.$$

U unutrašnjem čvoru $ij$ diferencijalna jednadžba postaje

$$\frac{u_{i\,j+1}-u_{i\,j}}{\tau} = c^2\,\frac{u_{i-1\,j} -2\,u_{i\,j} +u_{i+1\,j}}{h^2}.$$

Pomnožimo jednadžbu s $\tau,$ i stavimo

$$\sigma = \frac{\tau}{h^2}.$$

Tada jednadžba u čvoru $ij$ glasi

$$u_{i\,j+1} = c^2\,\sigma\,u_{i-1\,j} + (1 -2\,c^2\,\sigma)\,u_{i\,j} + c^2\,\sigma\,u_{i+1\,j}$$ (3.33)

Sada, za razliku od postupka kod jednadžbe ravnoteže, rješenje tražimo sukcesivno. Zahvaljujući upotrebljenoj diferencijskoj shemi,

najprije računamo $u_{i\,1}$ po formuli (3.33)

$$u_{i\,1} = c^2\,\sigma\,u_{i-1\,0} + (1 -2\,c^2\,\sigma)\,u_{i\,0} + c^2\,\sigma\,u_{i+1\,0},\qquad i=1,2,\ldots{},n-1.$$

Zatim formulom (3.33) za $j=2$ računamo $u_{i\,2},$

$$u_{i\,2} = c^2\,\sigma\,u_{i-1\,1} + (1 -2\,c^2\,\sigma)\,u_{i\,1} + c^2\,\sigma\,u_{i+1\,1},\qquad i=1,2,\ldots{},n-1$$

itd.

Da bi ovako definiran postupak bio korektan, on treba biti

konvergentan, tj. sa smanjivanjem koraka mora približno rješenje težiti k točnom rješenju, i

stabilan, tj. male promjene ulaznih podataka (greške zaokruživanja) ne smiju prouzročiti velike razlike u rješenju.

Može se pokazati da su oba uvjeta ispunjena, ako je

$$c^2\,\sigma\leqslant{}\frac{1}{2},$$   tj.$$\quad \tau\leqslant{}\frac{1}{2\,c^2}\,h^2.$$

Ovdje opisan postupak rješavanja, zasnovan na formuli (3.33) se zove eksplicitan. On je kao što se vidi vrlo jednostavan. Nedostatak postupka je u tome što se zbog potrebe stabilnosti i konvergencije mora korak po osi $t$ uzeti vrlo malen. Za $c=1$ i $h=0.1$ treba biti $\tau\leqslant{}0.005.$ Postoji implicitni postupak, kod kojeg nema zahtjeva na $\sigma.$ Međutim on vodi na rješavanje sustava jednadžbi, pa se tu pojavljuju drugi problemi, osobito kad su koraci maleni. Osim toga kod implicitnog postupka moramo unaprijed ograničiti vrijeme, kako bismo imali konačno mnogo čvorova, dok kod eksplicitnog postupka ne treba unaprijed ograničavati vrijeme.

Valna jednadžba

Problem koji sada rješavamo je

$$\frac{\textstyle{\partial^2 u(x,t)}}{\textstyle{\partial t^2}} = c^2\,\frac{\textstyle{\partial^2 u(x,t)}}{\textstyle{\partial x^2}},\quad x\in [0,\ell],\;t\in [0,\infty{}\rangle{},$$

$$u(0,t)=0,\quad u(\ell,t)=0,\quad t\in [0,\infty{}\rangle{},$$

$$u(x,0) = \alpha{}(x),\quad \frac{\textstyle{\partial u(x,0)}}{\textstyle{\partial t}} = \beta{}(x), \quad x\in [0,\ell].$$

Diskretizaciju područja učinimo kao kod jednadžbe provođenja. Razlika je u tome što u jednadžbi dolazi druga parcijalna derivacija po $t,$ i što imamo dodatni početni uvjet. Parcijalnu derivaciju iz početnog uvjeta možemo na razne načine aproksimirati diferencijama. Neka je na pr.

$$\frac{u_{i\,1}-u_{i\,0}}{\tau} = \beta_i.$$

U unutrašnjem čvoru $ij$ jednadžba je

$$\frac{u_{i\,j-1}-2\,u_{i\,j}+u_{i\,j+1}}{\tau^2} = c^2\,\frac{u_{i-1\,j} -2\,u_{i\,j} +u_{i+1\,j}}{h^2}.$$

Nakon množenja s $\tau^2,$ i stavljanja

$$\sigma = \frac{\tau}{h},$$

imamo eksplicitni postupak dan formulom

$$u_{i\,j+1} = -u_{i\,j-1} + c^2\,\sigma^2\,u_{i-1\,j} + 2\,(1-c^2\,\sigma^2)\,u_{i\,j} + c^2\,\sigma{}^2\,u_{i+1\,j}.$$ (3.34)

Vrijednosti u prva dva reda čvorova $(j=0,1)$ dobijemo iz početnih uvjeta. Vrijednosti u daljnjim redovima računamo sukcesivno pomoću formule (3.34). Za $\sigma{}\leqslant{}1$ postupak je stabilan.

Varijacijske metode

Ritzova metoda

Pogledajmo kako ona funkcionira kad se radi o rubnom problemu

$$- \Delta\,u = f,\quad \Omega$$

$$u\vert _{\partial\Omega} = 0.$$

Kao u jednodimenzionalnom slučaju izaberemo $n$ linearno nezavisnih funkcija $v_i, i=1,2,\ldots, n$ koje zadovoljavaju rubni uvjet. Rješenje se pretpostavi u obliku

$$u_n = \sum_{i=1}^n c_i\,v_i,$$

i neodređeni koeficijenti se odrede iz uvjeta da $u_n$ minimizira pripadni funkcional energije

$$F(u) = \frac{1}{2}\,\iint_{\Omega}\,({\rm grad\,}u)^2dxdy - \iint_{\Omega}\,f\,u\,dxdy.$$

Tako dobiveni $u_n$ leži u vektorskom prostoru razapetom s funkcijama $v_1,$ $v_2,$ $\ldots,$ $v_n.$ Rješenje ne mora ležati u tom prostoru, pa u tom slučaju $u_n$ nije točno već samo približno rješenje. No, što veći $n$ uzmemo, to je manja greška koju činimo prihvaćajući $u_n$ kao rješenje problema.

Prvi problem s kojim se susrećemo kod Ritzove metode je određivanje funkcija $v_i$ koje moraju zadovoljavati rubni uvjet. Ako područje nije dovoljno lijepo, mogu nastati problemi. Nakon što smo izabrali funkcije $v_i,$ pretpostavljeno rješenje uvrstimo u funkcional

$$F(u_n) = \frac{1}{2}\,\iint_{\Omega}... ...,}\sum_{i=1}^n c_i\,v_i)^2dxdy - \iint_{\Omega}\,f\,\sum_{i=1}^n c_i\,v_i\,dxdy$$

$$= \frac{1}{2}\,\iint_{\Omega}\, (\su... ...n c_i\,{\rm grad\,}v_i)^2dxdy - \sum_{i=1}^n c_i\,\iint_{\Omega}\, f\,v_i\,dxdy$$

$$= \frac{1}{2}\,\sum_{i=1}^n\sum_{j=1... ...v_i\cdot{\rm grad\,}v_j\,dxdy - \sum_{i=1}^n c_i\,\iint_{\Omega}\, f\,v_i\,dxdy$$

$$= \Phi(c_1,c_2,\ldots,c_n).$$

$\Phi$ je derivabilna funkcija od $n$ varijabli $c_1,c_2,\ldots,c_n,$ pa jednadžbe

$$\frac{\partial \Phi}{\partial c_j} =... ..., {\rm grad\,}v_i\cdot{\rm grad\,}v_j\,dxdy - \iint_{\Omega}\,f\,v_i\,dxdy = 0,$$

za $j=1,2,\ldots,n,$ predstavljaju nužan uvjet za ekstrem funkcije $\Phi$ u točki $(c_1,c_2,\ldots,c_n).$ Ovo je sustav od $n$ linearnih algebarskih jednadžbi od $n$ nepoznanica. Stavimo

$$K_{ij} = \iint_{\Omega}\,{\rm grad\,}v_i\cdot{\rm grad\,}v_j\,dxdy,\hspace{1cm}b_j = \iint_{\Omega}\, f\,v_j\,dxdy.$$

Sada se sustav može kratko zapisati

$$\boldsymbol{K}\,\boldsymbol{c} = \boldsymbol{b},$$

gdje je $\boldsymbol{K}=[K_{ij}],\boldsymbol{c}=[c_j],\boldsymbol{b}=[b_i].$

Nedostaci ove metode su u tome što je za proizvoljno područje teško naći funkcije $v_i,$ i u tome što je matrica $\boldsymbol{K},$ koja se inače zove matrica krutosti, puna matrica, tj. općenito je svaki njezin element različit od nule.

Metoda konačnih elemenata

Navedene nedostatke donekle ispravlja metoda konačnih elemenata. Područje $\Omega$ podijelimo jednostavnim likovima (trokutima četverokutima i sl.) na dijelove, elemente. Vrhove elemenata zovemo čvorovima. Podjela se vrši tako da niti jedan čvor ne leži na stranici nekog drugog čvora. Numeriramo elemente i čvorove. Čvorove treba pažljivo numerirati, da matrica sustava koji ćemo konačno dobiti bude što uža, tj. da pojas oko glavne dijagonale izvan kojeg su same nule bude što uži. Zatim odaberemo koordinatne funkcije $v_i, i=1,2,\ldots, n$ tako da $v_i$ u čvoru $i$ ima vrijednost $1,$ a u ostalim čvorovima vrijednost $0.$ Na elementima, koji nemaju čvor $i$ kao vrh, $v_i$ ima vrijednost $0.$ Na elementima, koji imaju čvor $i$ kao svoj vrh, funkciju $v_i$ definiramo kao polinom prvog stupnja ili viših stupnjeva ako to zahtijeva problem koji rješavamo. Objasnimo sada na jednom jednostavnom primjeru metodu konačnih elemenata za problem ravnoteže membrane.

Primjer 3.22   Neka na napetu kvadratnu membranu duljine stranice $1,$ učvršćene na rubu u ravnini $xy,$ djeluje vanjska sila u smjeru osi $z$ po iznosu jednaka $1.$ Treba naći ravnotežni položaj membrane.

Rješenje. Rubni problem za ovaj zadatak glasi, uz dodatnu pretpostavku da je $c=1,$

$$- \Delta\,u = 1,\quad \Omega$$

$$u = 0,\quad \partial{}\Omega{},$$

gdje je $\Omega$ kvadrat stranice $1,$ koji možemo uzeti u ravnini $xy$ tako da su mu vrhovi točke $(0,0),(1,0),(1,1),(0,1).$ Podijelimo na devet jednakih kvadrata i zatim svaki od tih kvadrata na dva trokuta. Time je $\Omega$ podijeljen na 18 elemenata (trokuta). Također imamo 16 čvorova, od toga 12 na rubu. To su točke

$$(0,0), (1/3,0), (2/3,0), (1,0), (0,1/3), (1/3,1/3), (2/3,1/3), (1,1/3),$$

$$(0,2/3), (1/3,2/3), (2/3,2/3), (1,2/3), (0,1), (1/3,1), (2/3,1), (1,1).$$

Na ovoj slici su brojevima označeni čvorovi, a brojevima u kružnicama konačni elementi (u ovom slučaju trokuti).

Funkcional energije u ovom slučaju je

$$F(w) = \frac{1}{2}\,\iint_{\Omega}\,({\rm grad\,}w)^2dxdy - \iint_{\Omega}\,w\,dxdy.$$

Za svaki čvor definiramo koordinatnu funkciju

$$v_1,v_2,v_3,\ldots{},v_{16},$$

tako da stavimo $v_i = 1$ u $i$-tom čvoru, a u ostalim čvorovima $v_i = 0.$ Na elementu $e$ koji ima $i$-ti čvor kao vrh stavimo

$$v_i^e(x,y) = \alpha{}_i^e + \beta{}_i^e\,x + \gamma{}_i^e\,y.$$

Za $i=6$ imamo $v_6^1,v_6^2,v_6^4,v_6^7,v_6^9,v_6^{10}.$ Nađimo na pr. $v_6^9.$ Zbog svojstva koordinatnih funkcija imamo

$$\alpha{}_6^9 + \beta{}_6^9\,\frac{1}{3} + \gamma{}_6^9\,\frac{1}{3} = 1$$

$$\alpha{}_6^9 + \beta{}_6^9\,\frac{2}{3} + \gamma{}_6^9\,\frac{1}{3} = 0$$

$$\alpha{}_6^9 + \beta{}_6^9\,\frac{2}{3} + \gamma{}_6^9\,\frac{2}{3} = 0.$$

Rješenje je $\alpha{}_6^9=2,\beta{}_6^9=-3,\gamma{}_6^9=0,$ pa je

$$v_6^9 = 2 - 3\,x.$$

Na isti način možemo naći da je

$$v_6^1 = 3\,y,\;v_6^2 = 3\,x,\;v_6^4 = 1 - 3\,x + 3\,y,\;v_6^7 = 1 + 3\,x - 3\,y,\;v_6^{10} = 2 - 3\,y,$$

a također i ostale koordinatne funkcije. Primijetimo da je zbog geometrijskih razloga

$$v_6^1 = v_7^3,\;v_6^2 = v_{10}^8,\;v_{10}^7 = v_{11}^9,\;v_7^4 = v_{11}^{10},\;v_6^7 = v_{11}^{15},\;v_6^4 = v_{11}^{12},$$

$$v_6^{10} = v_7^{12},\;v_6^9 = v_{10}^{15},\;v_7^{11} = v_{11}^{17},\;v_{10}^{16} = v_{11}^{18}.$$

Rješenje tražimo u obliku

$$u(x,y) = c_1\,v_1(x,y) + c_2\,v_2(x,y) + \cdots{} + c_{16}\,v_{16}(x,y).$$

U čvoru na rubu, na pr. prvom, $(0,0),$ funkcija $u$ se poništava. Dakle

$$0 = c_1\,v_1(0,0) + c_2\,v_2(0,0) + \cdots{} + c_{16}\,v_{16}(0,0).$$

U prvom čvoru funkcija $v_1$ prima vrijednost $1,$ a ostale funkcije primaju vrijednost $0.$ Tako imamo

$$0 = c_1.$$

Za svaki čvor na rubu možemo na taj način dobiti da je, zbog uvjeta na rubu, pripadni koeficijent jednak $0.$ Tako imamo

$$c_1 = c_2 = c_3 = c_4 = c_5 = c_8 = c_9 = c_{12} = c_{13} = c_{14} = c_{15} = c_{16} = 0.$$

Preostaje dakle

$$u(x,y) = c_6\,v_6(x,y) + c_7\,v_7(x,y) + c_{10}\,v_{10}(x,y) + c_{11}\,v_{11}(x,y).$$

Dakle zadatak je odrediti koeficijente $c_6,c_7,c_{10},c_{11}.$ To ćemo učiniti tako da u Bernoullijev princip (2.12) uvrstimo ovaj $u,$ i redom $v=v_i,$ gdje je $i=6,7,10,11.$ Tako imamo

$$\iint_{\Omega}\,{\rm grad\,}v_i\cdot{\rm grad\,}\sum_j c_j\,v_j\,dxdy = \iint_{\Omega}\, v_i\,dxdy,\quad i=6,7,10,11,$$

$$\sum_j c_j\,\iint_{\Omega}\,{\rm grad\,}v_i\cdot{\rm grad\,}v_j\,dxdy = \iint_{\Omega}\, v_i\,dxdy,\quad i=6,7,10,11.$$

Da bi se dobila jednadžba za čvor $i=6,$ treba izračunati

$$\iint_{\Omega}\,{\rm grad\,}v_6\cdot{\rm grad\,}v_j\,dxdy,\quad j=6,7,10,11,$$   i $$\quad \iint_{\Omega}\,v_6\,dxdy.$$

Integral po $\Omega$ je suma integrala po elementima. Zbog svojstava koordinatnih funkcija, integral po mnogim elementima iščezava. Na pr.

$$\iint_{\Omega}\,{\rm grad\,}v_6\cdot... ...rad\,}v_6\,dxdy = \iint_1 + \iint_2 + \iint_4 + \iint_7 + \iint_9 + \iint_{10}.$$

$${\rm grad\,}v_6^1 = \{0,3\},\;{\rm grad\,}v_6^2 = \{3,0\},\;{\rm grad\,}v_6^4 = \{-3,3\},$$

$${\rm grad\,}v_6^7 = \{3,-3\},\;{\rm grad\,}v_6^9 = \{-3,0\},\;{\rm grad\,} v_6^{10} = \{0,-3\}.$$

Tako je

$$\iint_1 = \iint_2 = \iint_9 = \iint_{10} = 9\,\frac{1}{18} = \frac{1}{2},\quad \iint_4 = \iint_7 = 18\,\frac{1}{18} = 1,$$

i prema tome

$$\iint_{\Omega}\,{\rm grad\,}v_6\cdot{\rm grad\,}v_6\,dxdy = 4.$$

Zatim

$$\iint_{\Omega}\,{\rm grad\,}v_6\cdot{\rm grad\,}v_7\,dxdy = \iint_4 + \iint_9.$$

$$v_7^4 = -1 + 3\,x,\quad v_7^9 = 3\,x - 3\,y,$$

$${\rm grad\,}v_7^4 = \{3,0\},\quad {\rm grad\,}v_7^9 = \{3,-3\},$$

pa je

$$\iint_4 = -9\,\frac{1}{18} = -\frac{1}{2},\quad \iint_9 = -9\,\frac{1}{18} = -\frac{1}{2},$$

i prema tome

$$\iint_{\Omega}\,{\rm grad\,}v_6\cdot{\rm grad\,}v_7\,dxdy = -1.$$

$$\iint_{\Omega}\,{\rm grad\,}v_6\cdot{\rm grad\,}v_{10}\,dxdy = \iint_7 + \iint_{10}.$$

$$v_{10}^7 = -1 + 3\,y,\quad v_{10}^{10} = - 3\,x + 3\,y,$$

$${\rm grad\,}v_{10}^7 = \{0,3\},\quad {\rm grad\,}v_{10}^{10} = \{-3,3\},$$

pa je

$$\iint_7 = -9\,\frac{1}{18} = -\frac{1}{2},\quad \iint_{10} = -9\,\frac{1}{18} = -\frac{1}{2},$$

i prema tome

$$\iint_{\Omega}\,{\rm grad\,}v_6\cdot{\rm grad\,}v_{10}\,dxdy = -1.$$

$$\iint_{\Omega}\,{\rm grad\,}v_6\cdot{\rm grad\,}v_{11}\,dxdy = \iint_9 + \iint_{10}.$$

$$v_{11}^9 = v_{10}^7 = -1 + 3\,y,\quad v_{11}^{10} = - 3\,x + 3\,y,$$

$${\rm grad\,}v_{11}^9 = \{0,3\},\quad {\rm grad\,}v_{11}^{10} = v_7^4 = \{3,0\},$$

pa je

$$\iint_9 = \iint_{10} = 0,$$

i prema tome

$$\iint_{\Omega}\,{\rm grad\,}v_6\cdot{\rm grad\,}v_{11}\,dxdy = 0.$$

Osim toga je

$$\iint_{\Omega}\,v_6\,dxdy = \iint_1 v_6^1 + \iint_2 v_6^2 + \iint_4 v_6^4 + \iint_7 v_6^7 + \iint_9 v_6^9 + \iint_{10} v_6^{10} = \frac{2}{27}.$$

Isti rezultat dobivamo i ako umjesto $v_6$ stavimo $v_7,v_{10}$ ili $v_{11}.$ Tako je jednadžba za čvor $i=6$

$$4\,c_6 - c_7 - c_{10} = \frac{2}{27}.$$

Također je očito

$$\iint_{\Omega}\,{\rm grad\,}v_i\cdot... ...iint_{\Omega}\,{\rm grad\,}v_6\cdot{\rm grad\,} v_6\,dxdy = 4,\qquad i=7,10,11.$$

$$\iint_{\Omega}\,{\rm grad\,}v_7\cdot... ...11}\,dxdy = \iint_{\Omega}\,{\rm grad\,}v_6\cdot{\rm grad\,} v_{10}\,dxdy = -1,$$

pa jednadžba za čvor $i=7$ glasi

$$- c_6 + 4\,c_7 - c_{11} = \frac{2}{27}.$$

$$\iint_{\Omega}\,{\rm grad\,}v_{10}\c... ...v_{11}\,dxdy = \iint_{\Omega}\,{\rm grad\,} v_6\cdot{\rm grad\,}v_7\,dxdy = -1,$$

pa jednadžbe za čvorove $i=10,11$ glase

$$- c_6 + 4\,c_{10} - c_{11} = \frac{2}{27}.$$

$$- c_7 - c_{10} + 4\,c_{11} = \frac{2}{27}.$$

Rješenje sustava je

$$c_6 = c_7 = c_{10} = c_{11} = \frac{1}{27}.$$

To su ujedno vrijednosti funkcije $u$ u čvorovima $i=6,7,10,11.$

Još o greškama

Cijeli brojevi

Za zapisivanje brojeva općenito, pa i u kompjuteru, koristi se niz znamenaka u nekom brojevnom sustavu. U tu svrhu treba odrediti bazu $\beta.$ Na pr. u dekadskom sustavu je $\beta=10,$ a u binarnom $\beta=2.$ Računati ne možemo s beskonačno mnogo znamenaka. To ne može niti kompjuter. Dapače, on je ograničen svojim fizičkim mogućnostima, hardwareom. Ograničenje se odnosi na broj znamenaka $p$ koje može upotrebiti za zapis jednog broja. Broj $p$ se zove preciznost. Uređeni par $(\beta,p)$ zove se format.

Cijeli brojevi se u kompjuteru zapisuju u formatu $(2,16)$ ili $(2,32).$ Dakle radi se o brojevima koji se zapisuju pomoću znamenki 0 i $1$ razmještenih na $16$ ili $32$ mjesta. Pri tom se prvo mjesto koristi za predznak. Broj 0 označava pozitivan broj, $1$ negativan. Pojedino mjesto na koje možemo upisati jednu znamenku zove se bit.

Imamo sljedeći algoritam za pretvaranje broja zapisanog u dekadskom sustavu u binarni.

Algoritam 13   Neka je $N$ dekadski zapis broja. Računamo

$$N_0 = N,\hspace{1cm}N_{k+1} = \frac{N_k-a_k}{2}$$    za $$k=0,1,2,\ldots{},$$

gdje je

$$a_k = \left\{ \begin{array}{ll} 1, &... ...{ako je $N_k$ neparan}, \\ 0, & \text{ako je $N_k$ paran}. \end{array}\right.$$

Račun provodimo sve dok ne dobijemo $N_k=0.$

Primjer 3.23   Treba za broj koji ima dekadski zapis $1998$ naći binarni zapis.

Rješenje.

$k$	$$0$$	$1$	$2$	$3$	$4$	$5$	$6$	$7$	$8$	$9$	$10$
$N_k$	$1998$	$999$	$499$	$249$	$124$	$62$	$31$	$15$	$7$	$3$	$1$
$a_k$	$$0$$	$1$	$1$	$1$	$$0$$	$$0$$	$1$	$1$	$1$	$1$	$1$

Dakle binarni zapis je $11111001110.$

Tako se u formatu $(2,16)$ dekadski broj $1998$ zapisuje kao

$$0000011111001110.$$

Segment određen donjom i gornjom granicom zapisa zove se rang. Budući da u formatu $(2,16)$ imamo $15$ bitova za upis znamenaka (jedan bit je rezerviran za predznak), i budući da se mogu koristiti samo znamenke 0 i $1,$ najveći broj koji tako možemo zapisati je

$$111111111111111 = 2^{14} + 2^{13} + \cdots + 2^1 + 2^0 = \frac{2^{15}-1}{2-1} = 2^{15}-1 = 32767.$$

Najmanji cijeli broj koji se tako može zapisati je $-32768.$ Tako je rang u ovom slučaju

$$[-32768,32767].$$

Ovakav zapis brojeva se zove zapis pomoću fiksnog zareza (fixed-point). U aritmetici foksnog zareza postoje greške, no mi ćemo se ovdje ograničiti samo na neke aspekte pojave grešaka kod zapisa i raračunanja s realnim brojevima.

Realni brojevi

Za kompjuterski zapis realnog broja obično se koristi format $(2,32)$ ili $(2,64).$ Za prvi format se kaže da se radi o jednostrukoj, a za drugi o dvostrukoj preciznosti.

U formatu $(\beta,p)$ realni broj se zapisuje u obliku

$$\pm d_0.d_{-1}d_{-2}d_{-3}\cdots d_{p} \times \beta^e,$$

gdje su $d$ znamenke $0\leqslant d_i<\beta,$ a $e$ je eksponent. Tako je u formatu $(2,32)$ $\beta=2,$ $p=32,$ a znamenke $d_i$ su 0 ili $1.$ U tom formatu prvi bit određuje predznak, a zadnjih osam eksponent, od kojih prvi određuje predznak eksponenta.

$$\pm d_0.d_{-1}d_{-2}d_{-3}\cdots d_{p} \times \beta^e,$$

Realan broj se uvijek zapisuje u normaliziranom obliku, tj. u obliku u kojem je $d_0\neq 0.$ To je zato da bi zapis broja bio jedinstven. Rang u ovom slučaju određuju minimalni i maksimalni eksponent. Maksimalni je

$$e_{max} = 2^0 + 2^1 + 2^2 + \cdots + 2^6 = 2^7-1 = 127,$$

dok je $e_{min}=-128.$ Rang određujemo tako da nađemo onaj $n$ za koji je

$$10^n = 2^{127}.$$

Ako logaritmiramo, imamo

$$n = 127\,\log 2 \approx 127\times 0.30103 \approx 38.$$

Dakle rang je

$$\left[10^{-38},10^{38}\right].$$

Ovakav zapis broja zovemo zapisom pomoću pomičnog zareza (floating-point).

Greške nastaju na dva načina. Jedan način da nastane greška je taj da se slijedom operacija dođe do broja koji nije u rangu. U tom slučaju kažemo da se dogodio overflow. O tome nećemo govoriti. Drugi način kako može doći do greške jeste kada želimo upisati ili kad operacijama dođemo do broja koji se ne može točno zapisati u formatu u kojem radi kompjuter. Na pr. dekadski broj $0.1$ se ne može točno zapisati u formatu s bazom $2$ bez obzira kolika bila preciznost.

Pretvorba dekadskog zapisa u binarni u slučaju realnog broja $0<z\leqslant 1$ dan je ovim algoritmom.

Algoritam 14   Neka je $z$ takav da je $0<z\leqslant{}1.$ Računamo nizove brojeva $(z_k)$ i $(a_{-k})$ iz formula

$$z_1 = z$$

$$% latex2html id marker 41886 a_{-k} = \left\{ \begin{array}{... ...\ 0, & \text{ ako je $2\,z_k\leqslant{}1,$} \end{array}\right.$$

$$z_{k+1} = 2\,z_k - a_{-k},\hspace{1cm}k=1,2,\ldots{}\ .$$

Tada je $z=0.a_{-1}a_{-2}a_{-3}\ldots{}$ binarni zapis broja $z.$

Pomoću ovog algoritma možemo utvrditi da binarni zapis dekadskog broja $0.1$ glasi

$$0.0001100110011001100\ldots\ .$$

Relativna greška

$$\frac{\vert z-z'\vert}{z}.$$

U zapisu pomoću pomičnog zareza je

$$z' = d.ddd\ldots{}d\times \beta{}^e,$$

pa je relativna greška

$$\frac{\vert z-d.ddd\ldots{}d\times \beta{}^e\vert}{z}.$$

Ako je $z'$ najbolja moguća aproksimacija broja $z,$ onda je apsolutna greška (brojnik ovog razlomka) najviše

$$\underbrace{0.000\ldots{}0}_{p\,\text{nula}}\frac{\beta{}}{2} \times \beta{}^e = \frac{\beta{}}{2}\,\beta{}^{-p} \times \beta{}^e.$$

Brojevi oblika $d.ddd\ldots{}d\times \beta{}^e$ se nalaze između $\beta{}^e$ i $\beta{}\,\beta{}^e.$ Prema tome najveća moguća relativna greška je

$$\frac{\frac{\beta{}}{2}\,\beta{}^{-p} \times \beta{}^e}{\beta{}^e} = \frac{\beta{}}{2}\,\beta{}^{-p} = \varepsilon{}.$$

Ovaj broj se zove kompjuterski epsilon. Najmanja moguća relativna greška je

$$\frac{\frac{\beta{}}{2}\,\beta{}^{-p} \times \beta{}^e}{\beta{}\,\beta{}^e} = \frac{1}{2}\,\beta{}^{-p}.$$

Tako je

$$\frac{1}{2}\,\beta{}^{-p} \leqslant{} \frac{\vert z-d.ddd\ldots{}d\times \beta{}^e\vert}{z} \leqslant{} \frac{\beta{}}{2}\,\beta{}^{-p}.$$

Osim relativnom greškom, grešku možemo mjeriti i jedinicom ulp (units in the last place). Ako je format $(\beta{},p),$ onda se greška u ulpima izražava formulom

$$\left\vert\frac{z}{\beta{}^e} - d.ddd\ldots{}d\right\vert\,\beta{}^{p-1}.$$

Primjer 3.24   Neka je format $(10,3),$ $z=0.0314,$ i $z'=3.12\times{}10^{-2}.$ Naći grešku u ulpima.

Rješenje. $d.dd = 3.12.$ Zatim $e=-2.$ Odatle

$$\left\vert\frac{0.0314}{10^{-2}} - 3.12\right\vert\,10^{2} = 0.02 \times{} 10^2 = 2 ulpa.$$

Greška u ulpima omogućava da odredimo broj nepouzdanih znamenki u približnoj vrijednosti $z'$ broja $z.$ Ako je greška $n$ ulpa, onada je broj nepouzdanih znamenki

$$\log_\beta{}n.$$

Primjer 3.25   U formatu $(10,3)$ treba izračunati $x-y,$ ako je $x=10.1,$ $y=9.93,$ i odrediti broj nepouzdanih znamenki približne vrijednosti.

Rješenje. Prilikom zbrajanja ili odbijanja manji broj svodimo na potenciju baze većega, i zatim zbrajamo samo signifikande

$$x = 10.1 = 1.01\times{}10^1$$   i$$\quad y = 9.93 = 0.993\times{}10^1.$$

Odatle

$$x' = 1.01\times{}10^1$$   i$$\quad y = 0.99\times{}10^1.$$

Tako je

$$x' - y' = 0.2\times{}10^1 = 2.00\times{}10^{-1},$$

dok je točna vrijednost

$$x - y = 0.17.$$

Greška u ulpima je prema tome

$$\left\vert\frac{0.17}{10^{-1}} - 2.00\right\vert\,10^2 = 0.3\times{}10^2 = 30$$ ulpa.

Budući da je

$$1 < \log 30 < 2,$$

nepouzdane su dvije znamenke.

U ovom primjeru je $x'=x,$ a greška zaokruživanja broja $y'$ je

$$\left\vert\frac{9.93}{10^0} - 9.90\right\vert\,10^2 = 3$$ ulpa.

Dakle u

$$y' = 9.90\times{}10^0$$

je jedna znamenka nepouzdana. Ipak prilikom odbijanja se broj nepouzdanih znamenki povećao. Dapače, kad se odbijaju bliski brojevi može se dogoditi da niti jedna znamenka u rezultatu nije pouzdana.

Taj problem se može ublažiti tako da se račun izvede s dodatnom znamenkom (guard digit).

Primjer 3.26   Riješimo primjer 3.25 pomoću dodatne znamenke.

Rješenje. Tada je

$$x' = 1.010\times{}10^1$$   i$$\quad y' = 0.993\times{}10^1,$$

$$x' - y' = 0.017\times{}10^1,$$

pa je

$$(x' - y')' = 0.017\times{}10^1 = 1.70\times{}10^{-1},$$

što je točna vrijednost.

Primjer 3.27   Neka je format $(10,3),$ $b=3.34,$ $a=1.22,$ i $c=2.28.$ Treba izračunati $b^2-4ac.$

Rješenje.

$$b^2 = 11.1556,\quad 4\,a\,c = 11.1264,\quad b^2 - 4\,a\,c = 0.0292,$$

$${b^2}' = 1.12\times{}10^1,\quad (4\,a\,c)' = 1.11\times{}10^1,\quad (b^2)' - (4\,a\,c)' = 0.1 = 1.00\times{}10^{-1}.$$

Dakle greška je

$$\left\vert 1.00 - \frac{0.0292}{10^{-1}}\right\vert\,10^2 = 70.8$$ ulpa.

Prema tome dvije znamenke su nepouzdane.

Ako se radi s dodatnom znamenkom, greška postaje samo $8$ ulpa.

Osim računanjem s dodatnom znamenkom, greška se može umanjiti rearanžiranjem formule.

Primjer 3.28   Ako je $b^2\gg 4ac,$ i $b>0,$ onda je $\sqrt{b^2-4ac}\approx b,$ pa formula

$$\frac{-b + \sqrt{b^2-4\,a\,c}}{2\,a}$$

sadrži odbijanje bliskih brojeva, što može dovesti do velike pogreške. U ovom slučaju možemo racionalizirati brojnik

$$\frac{-b + \sqrt{b^2-4\,a\,c}}{2\,a} = \frac{-b + \sqrt{b^2-4\,a\... ...t{b^2-4\,a\,c}}{-b - \sqrt{b^2-4\,a\,c}} = \frac{2\,c}{-b - \sqrt{b^2-4\,a\,c}}$$

Primjer 3.29   Heronova formula za računanje površine trokuta glasi

$$P = \sqrt{s(s-a)(s-b)(s-c)},$$

gdje je

$$s = \frac{a+b+c}{2}.$$

Neka je dan trokut takav da je $a\approx b+c.$ Tada ova formula daje loš rezultat za površinu. Znatno bolja je formula

$$P_1 = \frac{\sqrt{(a+(b+c))(c-(a-b))(c+(a-b))(a+(b-c))}}{4}$$

Primjer 3.30   Format je $(10,3).$ Treba izračunati $x^2-y^2,$ ako je $x=3.34,$ $y=3.33.$

Rješenje. Točan rezultat je

$$x^2 = 11.1556,\quad y^2 = 11.0889,\quad x^2-y^2 = 0.0667\ .$$

Računanje po formuli $x^2-y^2$ daje

$$(x^2)' = 1.12\times 10^1,\quad (y^2)' = 1.11\times 10^1,\quad (x^2)'-(y^2)' = 0.01\times 10^1 = 0.1,$$

$$((x^2)'-(y^2)')' = 1.00\times 10^{-1}.$$

Ovaj račun daje grešku od

$$\vert 1 - 0.667\vert\times 10^2 = 33$$ ulpa$$.$$

No, ako računamo po formuli $(x-y)(x+y),$ onda imamo

$$x-y = 0.01 \Rightarrow (x-y)' = 1.00\times 10^{-2},$$

$$x+y = 6.67 \Rightarrow (x+y)' = 6.67\times 10^{0}.$$

Tako je

$$(x-y)'\,(x+y)' = 6.67\times 10^{-2},$$

što je točan rezultat.

Primjer 3.31   Format je $(10,3).$ Treba izračunati $x+y_1+y_2+\cdots+y_{10},$ gdje je $x=1.00,$ $y_1=y_2=\cdots=y_{10}=0.00493.$

Rješenje. Točan rezultat je

$$x+y_1+y_2+\cdots+y_{10} = 1.0493\ .$$

Ako računamo po formuli

$$((\cdots(x + y_1) + y_2) + \cdots ) + y_{10}),$$

onda imamo

$$x' = 1.00\times 10^0,\qquad y_1' = 0.00\times 10^0,$$

jer se prilikom zbrajanja ili odbijanja brojevi poravnavaju po eksponentu najvećega, i zatim zaokružuju. Tako imamo

$$x' + y_1' = 1.00\times 10^0 \Rightarrow (x' + y_1')' = 1.00\times 10^0.$$

Nakon toga, na isti način

$$((x' + y_1')' + y_2')' = 1.00\times 10^0,$$

itd. Na kraju imamo

$$((\cdots(x' + y_1')' + y_2')' + \cdots )' + y_{10}')' = 1.00\times 10^0.$$

Ako računamo tako da najprije zbrojimo male brojeve, pa zatim dodamo $x,$ onda imamo

$$y_1 = \cdots = y_{10} = 4.93\times 10^{-3},$$

pa je

$$y_1' = \cdots = y_{10}' = 4.93\times 10^{-3}.$$

Odatle

$$((\cdots(y_1' + y_2')' + y_3')' + \cdots )' + y_{10}')' = 49.3\times 10^{-3} = 4.93\times 10^{-2} = 0.0493\times 10^0.$$

Da se ovaj broj pribroji broju $x=1.00\times 10^0,$ mora ga se svesti na potenciju s eksponentom $0,$ zaokružiti ga i zatim dodati $x.$ To znači

$$((\cdots(y_1' + y_2')' + y_3')' + \cdots )' + y_{10}')' + x' = 0.05\times 10^0 + 1.00\times 10^0 = 1.05\times 10^0,$$

pa je i

$$((\cdots(y_1' + y_2')' + y_3')' + \cdots )' + y_{10}')' + x' = 1.05\times 10^0.$$

Točno zaokružene operacije

Ako se koristi dodatna znamenka (guard digit) prilikom zaokruživanja, pa se zatim rezultat na kraju zaokruži, greška naravno može biti veća nego ako se točan rezultat na kraju zaokruži.

Ako se operacije izvode točno, pa se konačni rezultat zaokruži na najbliži kompjuterski realni broj, onda kažemo da je zaokruživanje točno. Međutim, kompjuteru treba reći koji je to najbliži broj. Da problem bude jasniji, ostanimo kod baze $10,$ i preciznosti $3.$ Za broj $3.123,$ ili $3.128$ jasno je kako ćemo ga zaokružiti. Ako je znamenka koja se odbacuje manja od $5,$ onda se prethodna znamenka ne mijenja. Ako je znamenka koja se odbacuje veća od $5,$ onda se prethodna znamenka uveća za jedan. Postavlja se pitanje kako zaokružiti broj $3.125$? Jedan način je da se prethodna znamenka uveća za jedan. To može dovesti do velike greške.

Primjer 3.32   Prema rekurzivnoj formuli

$$x_{n+1} = (x_n - y) + y$$

treba mnogo puta računati rezultat.

Rješenje. Radi jednostavnosti uzmimo format $(10,3),$ $x=1.00,$ $y=-0.555=-5.55\times 10^{-1}.$ Imamo $x_0=1.00,$

$$x_1 = (1.00 + 0.555)' - 0.555 = 1.56 - 0.555 = 1.005 = 1.01,$$

$$x_2 = (1.01 + 0.555)' - 0.555 = 1.57 - 0.555 = 1.015 = 1.02,$$

$$x_3 = (1.02 + 0.555)' - 0.555 = 1.58 - 0.555 = 1.025 = 1.03,$$

$$\vdots$$

Dobivamo sve veće brojeve. Tako na primjer

$$x_{845} = 9.45\ .$$

Drugi način bi bio da se brojevi, kod kojih je broj $5$ prva znamenka koja se odbacuje, u pola slučajeva zaokružuju odbacivanjem, a u drugih pola uvećavanjem prethodne znamenke za jedan. Pravilo može biti na pr. sljedeće.

-

Ako je znamenka ispred broja $5$ parna, onda se $5$ odbacuje.

-

Ako je znamenka ispred broja $5$ neparna, onda se $5$ odbacuje, a prethodna znamenka se uvećava za jedan.

Lako se vidi da primjenjujući ovo pravilo u prethodnom primjeru dobivamo točan rezultat.