Rješavanje običnih diferencijalnih jednadžbi

Problem početnog uvjeta (Cauchyjev problem)

Želimo riješiti diferencijalnu jednadžbu

$$y' = f(x,y),$$

na segmentu $[a,b]$ uz početni uvjet $y(a)=\alpha.$

Metode približnog rješavanja koje ćemo sada opisati osnivaju se na sljedećoj ideji. Podijelimo segment $[a,b]$ na $n$ jednakih dijelova

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

Duljina svakog podsegmenta je

$$h = \frac{b-a}{n}.$$

Brojeve $x_i=a+ih$ zovemo čvorovima, a broj $h$ zovemo korakom. Stavimo

$$y_i = y(x_i).$$

Cilj nam je odrediti $y_i$ za svaki $i=1,2,\ldots,n.$ To činimo tako da derivaciju zamijenimo odgovarajućom algebarskom aproksimacijom, kojom dolazimo do rekurzivne formule, pomoću koje računamo $y_{i+1}$ iz poznatog $y_i.$ Na taj način, rješenje, koje je neprekidna funkcija, zamjenjujemo konačnim brojem njezinih vrijednosti. Opisani postupak se zove diskretizacija. Očekujemo da će za dovoljno mali $h$ brojevi $y_i$ dovoljno dobro aproksimirati prave vrijednosti funkcije. Važno svojstvo koje diskretizacija treba imati jeste da s povećanjem $n$ dobivamo sve bolje aproksimacije, odnosno, preciznije, da brojevi $y_i$ teže prema pravim vrijednostima funkcije kada $n\rightarrow\infty.$

Eulerova metoda

Neka je dana diferencijalna jednadžba

$$y' = f(x,y),$$

na segmentu $[a,b]$ uz početni uvjet $y(a)=\alpha.$

Na temelju Taylorovog teorema srednje vrijednosti imamo

$$y_{i+1} = y(x_i+h) = y(x_i) + y'(x_i)\,h + \frac{1}{2}\,y''(x_i+t\,h)\,h^2.$$ (3.23)

Zanemarimo zadnji član, koji sadrži $h^2,$ pa dobivamo približnu vrijednost

$$y_{i+1} \approx y_i + y'_i\,h.$$

Odatle

$$y_{i+1} \approx y_i + h\,f(x_i,y_i).$$

Tako imamo sljedeći algoritam

Algoritam 11   (Eulerova metoda) Izaberemo dovoljno velik prirodni broj $n.$ Za zadani korak

$$h = \frac{b-a}{n}$$

računamo niz brojeva $(y_i)$ po formuli

$$y_{i+1} = y_i + h\,f(x_i,y_i),\hspace{1cm}i=0,1,2,\ldots,n-1,$$

s tim da je

$$x_i = a + i\,h,$$

a početna vrijednost $y_0$ je određena početnim uvjetom.

Mathematica program 6   (Eulerova metoda)

f[x_,y_]=; (* upisuje se f(x,y), ako je jednadzba y'=f(x,y) *)
pocuvj=;   (* pocetni uvjet u obliku {x0,y0} *)
kraj=;     (* drugi rub segmenta na kojem se trazi rjesenje *)
n=;        (* broj podsegmenata podjele *)
h=(kraj-pocuvj[[1]])/n;
x=pocuvj[[1]];
y=pocuvj[[2]];
rj={};
While[x-h<kraj,
rj=Append[rj,{x,y}]//N;
y=y+h f[x,y];
x=x+h];
rj

Primjer 3.17   Riješiti Eulerovom metodom Cauchyjev problem

$$y' = y,\qquad y(0) = 1$$

na segmentu $[0,1],$ dijeleći ga na deset dijelova.

Rješenje. Koristeći se ovim programom, dobivamo sljedeće točke grafa.

$$(0,1.),(0.1,1.1),(0.2,1.21),(0.3,1.331),(0.4,1.4641),(0.5,1.61051),$$

$$(0.6,1.77156),(0.7,1.94872),(0.8,2.14359),(0.9,2.35795),(1.,2.59374).$$

Rezultat se vidi na slici. Puna linija predstavlja graf točnog rješenja $e^x,$ a crtkana linija spaja nađene točke dijelovima pravaca.

Ova metoda je vrlo jednostavna, ali i vrlo gruba (ocjena pogreške je vrlo gruba) tako da se rijetko kada upotrebljava.

Metoda Runge-Kutta

U metodi Runge-Kutta se također računa $y_{i+1}$ pomoću već poznate vrijednosti $y_i.$ Taj račun je točniji, jer uzima u obzir i neke međuvrijednosti funkcije $f.$ Imamo sljedeći algoritam

Algoritam 12   (Metoda Runge-Kutta) Izaberemo dovoljno velik prirodni broj $n.$ Za zadani korak

$$h = \frac{b-a}{n}$$

računamo niz brojeva $(y_n)$ po formuli

$$y_{i+1} = y_{i} + \frac{1}{6}\,\left(k_1 + 2\,k_2 + 2\,k_3 + k_4\right),\qquad i=0,1,2,\ldots,n-1,$$
gdje je

$$k_1 = h\,f(x_i,y_i)$$

$$k_2 = h\,f\left(x_i + {\scriptstyle{\frac{h}{2}}},y_i + \scriptstyle{\frac{k_1}{2}}\right)$$

$$k_3 = h\,f\left(x_i + {\scriptstyle{\frac{h}{2}}},y_i + \scriptstyle{\frac{k_2}{2}}\right)$$

$$k_4 = h\,f(x_i + h,y_i + k_3).$$

Početna vrijednost $y_0$ je određena početnim uvjetom.

Radi bolje preglednosti prilikom računanja se može koristiti sljedeća tabela.

$i$	$x_i$	$y_i$	$h\,f$
$$0$$	$x_0=a$	$y_0$	$k_1$	$k_1$
	$x_0+\frac{h}{2}$	$y_0+\frac{k_1}{2}$	$k_2$	$2\,k_2$
	$x_0+\frac{h}{2}$	$y_0+\frac{k_2}{2}$	$k_3$	$2\,k_3$
	$x_0+h$	$y_0+k_3$	$k_4$	$k_4$
				$\frac{1}{6}\,\sum\approx{}y_1-y_0$
$1$	$x_1$	$y_1$	$k_1$	$k_1$
	$x_1+\frac{h}{2}$	$y_1+\frac{k_1}{2}$	$k_2$	$2\,k_2$
	$x_1+\frac{h}{2}$	$y_1+\frac{k_2}{2}$	$k_3$	$2\,k_3$
	$x_1+h$	$y_1+k_3$	$k_4$	$k_4$
				$\frac{1}{6}\,\sum\approx{}y_2-y_1$
$2$	$x_2$	$y_2$	$k_1$	$k_1$
	...	...	...	...

Mathematica program 7   (Metoda Runge-Kutta)

f[x_,y_]=;      (* zadavanje f(x,y), ako je y'=f(x,y) *)
x=;y=;          (* po\1etni uvjet *)
h=;             (* korak *)
rj={{x,y}};
br=0;           (* po\1etna vrijednost broja\1a *)
n=;             (* broj koraka *)
While[br<n,
k1=h f[x,y];
k2=h f[x+h/2,y+k1/2];
k3=h f[x+h/2,y+k2/2];
k4=h f[x+h,y+k3];
br=br+1;
x=x+h;
y=y+(k1+2 k2+2 k3+k4)/6;
rj=Append[rj,{x,y}]];
N[rj,10]        (* ispis rezultata s deset znamenaka *)

Primjer 3.18   Riješiti isti zadatak kao u primjeru 3.17 metodom Runge-Kutta.

Rješenje. Primijenimo li gornji program, koji slijedi algoritam metode Runge-Kutta, dobivamo sljedeće točke grafa rješenja.

$x$	$y$
$$0$$	$1.$
$0.1$	$1.105170833$
$0.2$	$1.221402571$
$0.3$	$1.349858497$
$0.4$	$1.49182424$
$0.5$	$1.648720639$
$0.6$	$1.822117962$
$0.7$	$2.013751627$
$0.8$	$2.225539563$
$0.9$	$2.459601414$
$1.$	$2.718279744$

Egzaktno rješenje je $y=e^x,$ ispisano s deset znamenaka

$x$	$y$
$$0$$	$1.$
$0.1$	$1.105170918$
$0.2$	$1.221402758$
$0.3$	$1.349858808$
$0.4$	$1.491824698$
$0.5$	$1.648721271$
$0.6$	$1.8221188$
$0.7$	$2.013752707$
$0.8$	$2.225540928$
$0.9$	$2.459603111$
$1.$	$2.718281828$

Odmah možemo uočiti mnogo bolju točnost nego kod metode Eulera.

Postoji cijela familija metoda zasnovanih na ideji da se podsegment podijeli na još manje dijelove kako bi se izračunala vrijednost u sljedećem čvoru. Sve one nose naziv Runge-Kutta metode. Opisani postupak se najčešće upotrebljava.

Rubni problem

Metoda konačnih diferencija

Neka je dan rubni problem

$$u''(x) - q(x)\,u(x) + f(x) = 0$$

$$u(0) = 0, \hspace{1cm}u'(\ell) = 0.$$

Podijelimo segment $[0,\ell{}]$ na $n$ jednakih podsegmenata.

$$0 = x_0 < x_1 < x_2 < \cdots < x_n = \ell{}.$$

Točke $x_0,x_1,x_2,\ldots,x_n$ se zovu čvorovi. Točke $x_1,x_2,\ldots,x_{n-1}$ se zovu unutrašnji čvorovi, a $x_0,x_n$ se zovu rubni čvorovi. Svaki segment ima duljinu $h=\frac{\ell{}}{n}.$ Broj $h$ se zove korak. On se izabire malen, da točnost bude veća, ali sa smanjivanjem koraka se javljaju drugi nepoželjni efekti, pa u odabiru koraka treba biti oprezan. Ideja metode se sastoji u tome da se u svakom čvoru diferencijalna jednadžba zamijeni odgovarajućom algebarskom jednadžbom i zatim riješi tako dobiveni sustav algebarskih jednadžbi. U tu svrhu treba derivacije zamijeniti odgovarajućim algebarskim aproksimacijama. Rješenje koje tako dobijemo predstavlja približne vrijednosti rješenja u čvorovima.

Uočimo jedan unutrašnji čvor $x_i.$ Koristit ćemo sljedeće oznake

$$u(x_i) = u_i,\hspace{1cm}u'(x_i) = u'_i,\hspace{1cm}u''(x_i) = u''_i.$$

Na isti način kao u 3.5.1 imamo

$$u_{i+1} \approx u_i + u'_i\,h.$$

Odatle

$$u'_i \approx \frac{u_{i+1}-u_i}{h}.$$

Ova aproksimacija derivacije se zove aproksimacija s desna.

Umjesto aproksimacije s desna ponekad se koriste aproksimacija s lijeva ili centralna aproksimacija. Polazeći od formule (3.23), u kojoj zamjenimo $h$ s $-h,$

$$u_{i-1} = u(x_i-h) = u_i - u'_i\,h + \frac{1}{2}\,u''(x_i-s\,h)\,h^2,$$

dobivamo, zanemarivanjem člana s $h^2,$ aproksimaciju s lijeva

$$u'_i \approx \frac{u_i-u_{i-1}}{h}.$$

Centralnu aproksimaciju dobivamo, ako Taylorove formule

$$u_{i+1} = u(x_i+h) = u_i + u'_i\,h + \frac{1}{2}\,u''_i\,h^2 + \frac{1}{3!}\,u'''(x_i+t\,h)\,h^3,$$

$$u_{i-1} = u(x_i-h) = u_i - u'_i\,h + \frac{1}{2}\,u''_i\,h^2 - \frac{1}{3!}\,u'''(x_i-s\,h)\,h^3.$$

za $h$ i $-h$ oduzmemo i podijelimo s $2\,h.$ Tada imamo

$$u'_i = \frac{u_{i+1}-u_{i-1}}{2\,h} + \frac{h^2}{3!}\,\left[u'''(x_i+t\,h) + u'''(x_i-s\,h)\right],$$

pa ako zanemarimo član s $h^2,$ dobivamo

$$u'_i \approx \frac{u_{i+1}-u_{i-1}}{2\,h}.$$

Za aproksimaciju druge derivacije, uzimamo

$$u_{i+1} = u(x_i+h) = u_i + u'_i\,h + \frac{1}{2}\,u''_i\,h^2 + \frac{1}{3!}\,u'''_i\,h^3 + \frac{1}{4!}\,u^{(iv)}(x_i+t\,h)\,h^4,$$

$$u_{i-1} = u(x_i-h) = u_i - u'_i\,h + \frac{1}{2}\,u''_i\,h^2 - \frac{1}{3!}\,u'''_i\,h^3 + \frac{1}{4!}\,u^{(iv)}(x_i-s\,h)\,h^4.$$

Zbrojimo ove jednakosti

$$u_{i+1} + u_{i-1} = 2\,u_i + u''_i\,h^2 + \frac{h^4}{4!}\,\left[u^{(iv)}(x_i+t\,h) + u^{(iv)}(x_i-s\,h)\right].$$

Podijelimo s $h^2,$ i izračunamo $u''_i$

$$u''_i = \frac{u_{i+1}-2\,u_i+u_{i-1}}{h^2} - \frac{h^2}{4!}\,\left[u^{(iv)}(x_i+t\,h) + u^{(iv)}(x_i-s\,h)\right].$$

Zanemarimo zadnji član s faktorom $h^2,$ i dobivamo

$$u''_i \approx{}\frac{u_{i+1}-2\,u_i+u_{i-1}}{h^2}.$$

U čvoru $x=x_i$ diferencijalna jednadžba iz rubnog problema glasi

$$u''_i - q(x_i)\,u_i + f(x_i) = 0.$$

Zamijenimo li drugu derivaciju s njezinom aproksimacijom, dobijemo algebarsku jednadžbu za $i$-ti čvor

$$\frac{u_{i+1}-2\,u_i+u_{i-1}}{h^2} - q_i\,u_i + f_i = 0.$$

Pomnožimo jednadžbu s $-h^2,$ pa imamo

$$-u_{i-1} + (2+h^2\,q_i)\,u_i - u_{i+1} = h^2\,f_i,\hspace{1cm} i=1,2,\ldots,n-1.$$

U rubnim čvorovima imamo zadan rubni uvjet. Tako je na lijevom rubu

$$u_0 = 0,$$ (3.24)
a na desnom

$$\frac{u_n-u_{n-1}}{h} = 0,$$
tj.

$$u_n = u_{n-1}.$$ (3.25)

Na taj način smo dobili sustav od $n-1$-ne linearne algebarske jednadžbe. Rubni uvjeti imaju utjecaja samo na prvu jednadžbu $(i=1),$ koja postaje, zbog (3.24),

$$(2+h^2\,q_1)\,u_1 - u_{2} = h^2\,f_1,$$

i na zadnju jednadžbu $(i=n-1),$ koja, zbog (3.25), postaje

$$-u_{n-2} + (1+h^2\,q_{n-1})\,u_{n-1} = h^2\,f_{n-1}.$$

Sustav možemo zapisati u matričnom obliku. Ako uzmemo prirodan poredak jednadžbi počevši od čvora $1$ preko $2$ sve do $n-1,$ imamo vektor nepoznanica

$$% latex2html id marker 40746 u = \left[ \begin{array}{cccc} u_1 & u_2 & \cdots{} & u_{n-1} \end{array}\right]^T,$$

matricu sustava

$$% latex2html id marker 40748 A(h) = \left[ \begin{array}{ccc... ...& 0 & 0 & \cdots{} & -1 & (1+h^2\,q_{n-1}) \end{array}\right],$$

i vektor desne strane

$$% latex2html id marker 40750 b(h) = \left[ \begin{array}{ccc... ...\,h^2 & f_2\,h^2 & \cdots{} & f_{n-1}\,h^2 \end{array}\right]^T$$

Tada sustav možemo zapisati matrično

$$A(h)\,u = b(h).$$

Za kvadratnu matricu $A=[a_{ij}]$ kažemo da je striktno dijagonalno dominantna ako je

$$\sum_{\substack{i=0\\ i\neq k}}^n\vert a_{i\,j}\vert< \vert a_{i\,i}\vert,\qquad i=1,2,\ldots,n.$$

Može se pokazati da iz striktne dijagonalne dominantnosti matrice $A$ slijedi njezina regularnost. Ako je $q(x)>0,\forall x,$ onda je matrica $A(h)$ striktno dijagonalno dominantna, pa je regularna. Prema tome postoji jedinstveno rješenje.

Primjer 3.19   Riješiti metodom konačnih diferencija sljedeći rubni problem

$$\ln (1 + x)\,y(x) + x\,y'(x) + y''(x) = 0,$$

$$y(0) = 1,\qquad y'(1) = 0.$$

Rješenje. Podijelimo segment $[0,1]$ na deset jednakih dijelova, tako da je korak $h=0.1.$ Čvorovi su $x_i=i\,h,\;i=0,1,2,\ldots,10.$ Stavimo $y(x_i)=y_i.$ Za prvu derivaciju upotrebimo desnu aproksimaciju, osim na desnom rubu, gdje uzmemo lijevu. U svakom unutrašnjem čvoru zamijenimo derivacije odgovarajućim aproksimacijama. Dobivamo sustav jednadžbi

$$y_0 = 1$$

$$100\,y_0 - 200.905\,y_1 + 101\,y_2 = 0$$

$$100\,y_1 - 201.818\,y_2 + 102\,y_3 = 0$$

$$100\,y_2 - 202.738\,y_3 + 103\,y_4 = 0$$

$$100\,y_3 - 203.664\,y_4 + 104\,y_5 = 0$$

$$100\,y_4 - 204.595\,y_5 + 105\,y_6 = 0$$

$$100\,y_5 - 205.53\,y_6 + 106\,y_7 = 0$$

$$100\,y_6 - 206.469\,y_7 + 107\,y_8 = 0$$

$$100\,y_7 - 207.412\,y_8 + 108\,y_9 = 0$$

$$100\,y_8 - 208.358\,y_9 + 109\,y_{10} = 0$$

$$-10\,y_9 + 10\,y_{10} = 0$$

čije rješenje je

$${{y_0} = {1.}},\quad {{y_1} = {1.05352}},\quad {{y_2} = {1.10551}},\quad {{y_3} = {1.1545}},$$

$${{y_4} = {1.19913}},\quad {{y_5} = {1.23816}},\quad {{y_6} = {1.27055}},\quad {{y_7} = {1.29548}},$$

$${{y_8} = {1.31235}},\quad {{y_9} = {1.32083}},\quad {{y_{10}} = {1.32083}}.$$

Točke rješenja na sljedećoj slici spojene su ravnim linijama.

Ritzova metoda

Iz varijacijskog principa slijedi da se umjesto rješavanja jednadžbe, rješenje rubnog problema može dobiti rješavanjem problema minimizacije odgovarajućeg funkcionala.

Ima raznih metoda koje koriste varijacijsku formulaciju. Jedna od njih je Ritzova. Pogledajmo kako ona funkcionira kad se radi o rubnom problemu

$$(p(x)\,u'(x))' +f(x) = 0,\hspace{1cm}u(0) = 0,\;\;u'(\ell{}) = 0,$$

gdje je $p(x)>0$ za svaki $x.$

Izaberemo $n$ linearno nezavisnih funkcija $v_i, i=1,2,\ldots, n$ koje zadovoljavaju Dirichletov homogeni rubni uvjet. Dakle $v_i(0)=0.$ Rješenje se pretpostavi u obliku

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

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

$$F(u) = \frac{1}{2}\,\int_0^{\ell}\,p\,{u'}^2\,dx - \int_0^{\ell}\,f\,u\,dx.$$

Tako dobiveni $u_n$ leži u vektorskom prostoru razapetom s funkcijama $v_1,v_2,\ldots,v_n,$ tj. u vektorskom prostoru svih linearnih kombinacija funkcija $v_1,v_2,\ldots,v_n.$ Rješenje ne mora biti linearna kombinacija tih funkcija, pa u tom slučaju $u_n$ nije točno već samo približno rješenje. No, grešku koju tom pretpostavkom činimo, možemo umanjiti uzimanjem većeg $n.$

Prvi problem s kojim se susrećemo kod Ritzove metode je određivanje funkcija $v_i,$ koje se često zovu koordinatne funkcije. One se obično biraju u skladu s problemom koji se rješava. Na pr. ako iz fizikalnih razloga očekujemo periodičko rješenje, onda ćemo takvima pretpostaviti i funkcije $v_i.$ Inače se za koordinatne funkcije mogu koristiti polinomi.

Nakon što smo izabrali funkcije $v_i,$ pretpostavljeno rješenje uvrstimo u funkcional

$$F(u_n) = \frac{1}{2}\,\int_0^{\ell}\,p\,\left(\sum_{i=1}^n c_i\,v'_i\right)^2\,dx - \int_0^{\ell}\,f\,\sum_{i=1}^n c_i\,v_i\,dx$$

$$= \frac{1}{2}\,\sum_{i=1}^n\sum_{j=1}^n c_i\,c_j\,\int_0^{\ell}\,p\,{v'}_i\, v'_j\,dx - \sum_{i=1}^n c_i\,\int_0^{\ell}\, f\,v_i\,dx$$

$$= \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_i} = \sum_{j=1}^n c_j\,\int_0^{\ell}\, p\,v'_i\, v'_j\,dx - \int_0^{\ell}\,f\,v_i\,dx = 0,$$

za $i=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 s $n$ nepoznanica. Ako stavimo

$$K_{i\,j} = \int_0^{\ell}\,p\,v'_i\,v'_j\,dx,\hspace{1cm}b_i = \int_0^{\ell}\, f\,v_i\,dx,$$

dobivamo sustav jednadžbi

$$\sum_{j=1}^n K_{i\,j}\,c_j = b_i,\qquad i=1,2,\ldots,n.$$

On se može matrično zapisati

$$K\,\boldsymbol{c} = \boldsymbol{b},$$ (3.26)

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

Matrica $K,$ koja se inače zove matrica krutosti, je očito simetrična, jer je

$$K_{i\,j} = \int_0^{\ell}\,p\,v'_i\,v'_j\,dx = \int_0^{\ell}\,p\,v'_j\,v'_i\,dx = K_{ji}.$$

Ona je i pozitivno definitna. Doista,

$$K\,\boldsymbol{c}\cdot{}\boldsymbol{c} = \left[\sum_{j=1}^n K_{ij}\,c_j\right]\cdot{}[c_i] = \sum_{i=1}^n \sum_{j=1}^n K_{ij}\,c_j\,c_i =$$

$$\sum_{i=1}^n \sum_{j=1}^n \int_0^{\ell}\, p\,v'_i\,v'_j\,dx\,c_j\... ...i\,c_j\,v'_j\,dx = \int_0^{\ell}\,p\,\left(\sum_{i=1}^n c_i\,v'_i\right)^2\,dx.$$

Budući da je $p(x)>0,$ za svaki $x,$ ovaj integral se može poništavati samo tako da bude

$$\left(\sum_{i=1}^n c_i\,v'_i(x)\right)^2 = 0,\quad \forall x\in [0,\ell],$$

odnosno

$$\sum_{i=1}^n c_i\,v'_i(x) = \left(\sum_{i=1}^n c_i\,v_i(x)\right)' = 0,\quad \forall x\in [0,\ell].$$

Odatle

$$\sum_{i=1}^n c_i\,v_i(x) = C,\quad \forall x\in [0,\ell],$$

a kako je $\sum_{i=1}^n c_i\,v_i(0) = 0,$ slijedi $C=0,$ odnosno

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

Zbog linearne nezavisnosti koordinatnih funkcija $v_i,$ slijedi da je $c_i=0$ za svaki $i,$ tj. $\boldsymbol{c}=\textbf{0}.$ Dakle $K\boldsymbol{c}\cdot\boldsymbol{c}>0,$ za $\boldsymbol{c}\neq\textbf{0},$ pa je matrica $K$ pozitivno definitna.

Kao što smo vidjeli, takva matrica se može dijagonalizirati i ima pozitivne vlastite vrijednosti. To znači da je regularna. Zaista, neka su njezine vlastite vrijednosti $\lambda{}_1,\lambda{}_2,\ldots,\lambda{}_n.$ Tada je

$$\det{}K = \det{}(X^{-1}[\lambda{}_i\,\delta_{ij}]\,X) = \det{}X^{... ...i\,\delta_{ij}]\,\det{}X = \lambda{}_1\,\lambda{}_2\,\cdots{}\,\lambda{}_n > 0.$$

Odatle slijedi da jednadžba (3.26) ima jedno i samo jedno rješenje.

Nedostak ove metode je u tome što je matrica $K$ puna matrica, tj. općenito je svaki njezin element različit od nule.

Primjer 3.20   Riješimo Ritzovom metodom sljedeći rubni problem

$$((1+x^2)\,u'(x))' + x(x-1)=0,$$

$$u(0)=0,\qquad u'(1)=0.$$

Rješenje. Za koordinatne funkcije uzmimo polinome. Budući da na lijevom rubu imamo homogen Dirichletov uvjet, polinomi se moraju poništavati u nuli. Neka su, dakle, koordinatne funkcije

$$v_i(x) = x^i,\qquad i=1,2,3,4.$$

Rješenje pretpostavljamo u obliku

$$u_4(x) = {c_1}\,{v_1}(x) + {c_2}\,{v_2}(x) + {c_3}\,{v_3}(x) + {c_4}\,{v_4}(x).$$

Kad s $u_4$ uđemo u funkcional, i njegove derivacije po koeficijentima $c_1,c_2,c_3,c_4$ izjednačimo s nulom, dobivamo sustav jednadžbi

$$1.33333\,{c_1} + 1.5\,{c_2} + 1.6\,{c_3} + 1.66667\,{c_4} = -0.0833333$$

$$1.5\,{c_1} + 2.13333\,{c_2} + 2.5\,{c_3} + 2.74286\,{c_4} = -0.05$$

$$1.6\,{c_1} + 2.5\,{c_2} + 3.08571\,{c_3} + 3.5\,{c_4} = -0.0333333$$

$$1.66667\,{c_1} + 2.74286\,{c_2} + 3.5\,{c_3} + 4.06349\,{c_4} = -0.0238095$$

Odavde dobijemo za koeficijente

$$c_1 = -0.171183,\quad c_2 = 0.038559,\quad c_3 = 0.141054,\quad c_5 = -0.0831691.$$

pa je rješenje

$$u_4(x) = -0.171183\,x + 0.038559\,{x^2} + 0.141054\,{x^3} - 0.0831691\,{x^4}.$$

Točno rješenje je

$${{u(x)} = {{\frac{x}{2}} - {\frac{{x^2}}{6}} - {\frac{2\,{\rm Arctg}\,x}{3}} + {\frac{\ln (1 + {x^2})}{6}}}}.$$

Greška računata u točkama $0,.1,.2,.3,.4,.5,.6,.7,.8,.9,1$ iznosi

$$0,-0.00014589,0.0000281161, 0.000191892,0.000195642,0.000056621,$$

$$-0.000107126,-0.000166607, -0.0000700444,0.0000853256, 2.45198\,{{10}^{-6}}.$$

Metoda konačnih elemenata

Metoda konačnih elemenata je modifikacija Ritzove metode u tom smislu da što više elemenata matrice krutosti bude jednako nuli. Budući da su elementi matrice krutosti

$$K_{i\,j} = \int_0^{\ell}\,p\,v'_i\,v'_j\,dx,$$

funkcije $v_i$ treba birati tako da međusobni produkti njihovih derivacija budu nulfunkcije u što više slučajeva.

Pogledajmo na primjeru kako se to radi. Neka je dan rubni problem

$$(p(x)\,u'(x))' +f(x) = 0,\qquad u(0) = 0,\;\;u'(\ell{}) = 0.$$

Podijelimo segment $[0,\ell{}]$ na $n$ jednakih dijelova. Time smo dobili čvorove

$$x_i = i\,h,\qquad i=0,1,2,\ldots,n,$$

gdje je

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

korak diskretizacije. Rješenje tražimo u obliku

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

gdje su $v_i$ koordinatne funkcije definirane formulom

$$v_i(x) = \begin{cases}\frac{x}{h} - ... ...o je }x_i\leqslant{}x\leqslant{}x_{i+1} \\ [1mm] 0,& \text{inače}, \end{cases}$$ (3.27)

za $i=1,2,\ldots,n-1.$ Koordinatne funkcije $v_0$ i $v_n$ su definirane ovako

$$v_0(x) = \begin{cases}-\frac{x}{h} +... ...ako je }x_{0}\leqslant{}x\leqslant{}x_1 \\ [1mm] 0,& \text{inače}, \end{cases}$$ (3.28)

$$v_n(x) = \begin{cases}\frac{x}{h} - ... ...o je }x_{n-1}\leqslant{}x\leqslant{}x_n \\ [1mm] 0,& \text{inače}, \end{cases}$$ (3.29)

Rješenje će biti oblika

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

Imamo

$$v_i(x_j) = \delta_{i\,j},$$

pa je prema tome

$$u(x_j) = \sum_{i=0}^{n} c_i\,v_i(x_j) = \sum_{i=0}^{n} c_i\,\delta_{i\,j} = c_j.$$ (3.30)

Dakle, neodređeni koeficijenti $c_i$ su upravo vrijednosti približnog rješenja u odgovarajućim čvorovima. Budući da je uvjet na lijevom rubu homogen, imamo

$$u(0) = u(x_0) = c_0 = 0,$$

pa će, prema tome, rješenje biti oblika

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

Budući da na desnom rubu imamo prirodni (Neumannov) uvjet, nema nikakvih daljnjih zahtjeva na oblik rješenja.

Iz formule (3.27) slijedi

$$v_i'(x) = \begin{cases} \frac{1}{h}... ... je }x_i\leqslant{}x\leqslant{}x_{i+1} \\ [1mm] 0,& \text{inače}. \end{cases}$$

Odatle dobivamo

$$K_{i\,j} = \int_0^{\ell}\,p\,v'_i\,v... ...x_{i+1}} p(x)\,dx,& \text{ako je }j=i+1 \\ [2mm] 0,& \text{inače}. \end{cases}$$ (3.31)

Dalje radimo isto kao kod Ritzove metode, da bismo na kraju dobili sustav linearnih algebarskih jednadžbi

$$\sum_{j=1}^n K_{i\,j}\,c_j = b_i,\qquad i=1,2,\ldots,n.$$

Elementi $K_{i\,j}$ za fiksni $i$ čine $i$-ti redak matrice krutosti $K.$ Formula (3.31) nam pokazuje da u jednom retku matrica $K$ ima najviše tri elementa različita od nule. To znatno pojednostavnjuje rješavanje sustava.

Primjer 3.21   Riješimo primjer 3.20 metodom konačnih elemenata.

Rješenje. Budući da je uvjet na desnom rubu homogen i prirodan, njega ne treba uzimati u obzir kod izbora koordinatnih funkcija. Uvjet na lijevom rubu nam kaže da treba isključiti funkciju $v_0.$ Dakle rješenje treba tražiti kao

$$u_{10} = {c_1}\,{v_1} + {c_2}\,{v_2} + {c_3}\,{v_3} + {c_4}\,{v_... ..._6}\,{v_6} + {c_7}\,{v_7} + {c_8}\,{v_8} + {c_9}\,{v_9} + {c_{10}}\,{v_{10}}.$$

Kad se izračunaju $K_{ij}$ i $b_i,$ dobije se sljedeći sustav jednadžbi

$$20.2667\,{c_1} - 10.2333\,{c_2} = 0.00883333,$$

$$-10.2333\,{c_1} + 20.8667\,{c_2} - 10.6333\,{c_3} = 0.0158333,$$

$$-10.6333\,{c_2} + 21.8667\,{c_3} - 11.2333\,{c_4} = 0.0208333,$$

$$-11.2333\,{c_3} + 23.2667\,{c_4} - 12.0333\,{c_5} = 0.0238333,$$

$$-12.0333\,{c_4} + 25.0667\,{c_5} - 13.0333\,{c_6} = 0.0248333,$$

$$-13.0333\,{c_5} + 27.2667\,{c_6} - 14.2333\,{c_7} = 0.0238333,$$

$$-14.2333\,{c_6} + 29.8667\,{c_7} - 15.6333\,{c_8} = 0.0208333,$$

$$-15.6333\,{c_7} + 32.8667\,{c_8} - 17.2333\,{c_9} = 0.0158333,$$

$$-17.2333\,{c_8} + 36.2667\,{c_9} - 19.0333\,{c_{10}} = 0.00883333,$$

$$-19.0333\,{c_9} + 19.0333\,{c_{10}} = 0.00158333$$

Iz ovog sustava se izračunaju koeficijenti

$${c_1} = 0.00883333,\;{c_2} = 0.0158333, {c_3} = 0.0208333,\;{c_4} = 0.0238333, {c_5} = 0.0248333,$$

$${c_6} = 0.0238333, {c_7} = 0.0208333,\;{c_8} = 0.0158333, {c_9} = 0.00883333,\;{c_{10}} = 0.00158333.$$

Budući da je rješenje između dva susjedna čvora linearna kombinacija polinoma prvog stupnja, rezultat će opet biti polinom prvog stupnja, tj. graf će biti dio pravca. Prema (3.30) rješenje je na rubu jednako odgovarajućim koeficijentima. Tako je za grafički prikaz rješenja dovoljno nanijeti točke $(x_i,c_i)$ kojima je dodana točka $(0,0)$ na lijevom rubu

$$(0,0),\,(0.1,0.00883333),\;(0.2,0.0158333),\;(0.3,0.0208333),\; (0.4,0.0238333),\;(0.5,0.0248333),$$

$$(0.6,0.0238333),\; (0.7,0.0208333),\;(0.8,0.0158333),\; (0.9,0.00883333),\;(1.,0.00158333),$$

i zatim te točke spojiti ravnim linijama.

Galerkinova metoda

Galerkinova metoda ili preciznije metoda Bubnova-Galerkina osniva se na sljedećoj jednostavnoj činjenici. Neka je dan vektorski prostor $V,$ baza $\boldsymbol{v}_1,\boldsymbol{v}_2,\ldots,\boldsymbol{v}_n,\ldots$ u $V$ i skalarni produkt u $V.$ Tada je $\boldsymbol{x}\in V$ nulvektor, ako i samo ako je $\boldsymbol{x}$ okomit na $\boldsymbol{v}_i$ za svaki $i.$ Ta tvrdnja je jasna, ako za $V$ uzmemo vektorski prostor radijvektora u prostoru ili ravnini.

Ono što nas zanima je rubni problem, na pr.

$$(p\,u')' + f = 0,\hspace{1cm}u(0)=a,\quad u'(\ell)=b.$$

Rješenje tražimo u skupu $D$ funkcija koje su klase $C^2([0,\ell]),$ i koje zadovoljavaju rubne uvjete. Taj skup nije vektorski prostor, jer linearne kombinacije takvih funkcija više ne zadovoljavaju ove rubne uvjete. Međutim, ako su rubni uvjeti homogeni, onda skup $D$ jeste vektorski prostor. Homogenizacijom rubnih uvjeta možemo svaki rubni problem svesti na problem s homogenim uvjetima.

Promatrajmo dakle rubni problem

$$(p\,u')' + f = 0,\hspace{1cm}u(0)=0,u'(\ell)=0.$$

U skupu $D,$ koji je sada vektorski prostor, definirajmo skalarni produkt

$$u\cdot{}v = \int_0^{\ell}\,u(x)\,v(x)\,dx.$$

Izaberimo u $D$ linearno nezavisne funkcije

$$\varphi_1,\varphi_2,\varphi_3,\ldots$$

tako da čine bazu u $D.$ U pravilu funkcija $v_i$ ima beskonačno mnogo. Neka je $u$ rješenje rubnog problema. Tada je $u\in D,$ i prema tome postoje brojevi $c_1,c_2,c_3,\ldots$ takvi da je

$$u(x) = \sum_{i=1}^{\infty} c_i\,v_i(x).$$

U diferencijalnoj jednadžbi rubnog problema, funkcija $(p\,u')' + f$ je jednaka nulfunkciji. To možemo na drugi način iskazati tako da zahtijevamo da je funkcija $(p\,u')' + f$ okomita na $v_i$ za svaki $i,$ pa nepoznati koeficijenti $c_i$ moraju zadovoljavati sljedeće jednadžbe

$$\int_0^{\ell}\,\left[\left(p(x)\,\left(\sum_{i=1}^{\infty} c_i\,v_i(x)\right)'\right)' + f(x)\right]\,v_j(x)\,dx = 0,\qquad j=1,2,3,\ldots\ .$$

Problem je u tome što je to beskonačno mnogo jednadžbi s beskonačno mnogo nepoznanica. Zato uzmemo konačno mnogo funkcija iz baze

$$v_1,v_2,\ldots,v_n,$$

i rješenje pretpostavimo u obliku

$$u_n = c_1\,v_1+c_2\,v_2+\ldots+c_n\,v_n = \sum_{j=1}^n c_j\,v_j,$$

Nepoznate koeficijente $c_1,c_2,\ldots,c_n$ određujemo iz sustava jednadžbi

$$\int_0^{\ell}\,\left[\left(p(x)\,\left(\sum_{j=1}^n c_j\,v_j(x)\right)'\right)' + f(x)\right]\,v_i(x)\,dx = 0,\qquad i=1,2,3,\ldots,n.$$

Ovaj sustav jednadžbi možemo prepisati u obliku

$$-\sum_{j=1}^n\,c_j\,\int_0^{\ell}\,\left(p(x)\,v'_j(x)\right)'\,v_i(x)\,dx = \int_0^{\ell}\,f(x)\,v_i(x)\,dx,\qquad i=1,2,3,\ldots,n.$$

Stavimo

$$K_{i\,j} = -\int_0^{\ell}\,\left(p(x)\,v'_j(x)\right)'\,v_i(x)\,dx,\qquad b_i = \int_0^{\ell}\,f(x)\,v_i(x)\,dx.$$

Tada sustav poprima oblik

$$\sum_{j=1}^n\,K_{i\,j}\,c_j = b_i,\qquad i=1,2,3,\ldots,n,$$

odnosno matrično

$$K\,\boldsymbol{c} = \boldsymbol{b},$$ (3.32)

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

U ovom slučaju se u formuli za $K_{ij}$ može jednom parcijalno integrirati pa, uzevši u obzir rubne uvjete, imamo

$$K_{i\,j} = - \left. p(x)\,v'_j(x)\,v_i(x)\right\vert _0^{\ell} + ... ...\ell}\, p(x)\,v'_j(x)\,v'_i(x)\,dx = \int_0^{\ell}\,p(x)\,v'_j(x)\,v'_i(x)\,dx.$$

Na taj način sustav jednadžbi (3.32) postaje identičan onome kod Ritzove metode. To se događa ako rubni problem ispunjava određene uvjete, o čemu ovdje nećemo detaljnije govoriti.

Istaknimo ovdje bitnu razliku u ideji između Ritzove i Galerkinove metode. Nužan uvjet za primjenu Ritzove metode je bila egzistencija varijacijske formulacije rubnog problema u kojoj se pojavljuje funkcional energije, dok za Galerkinovu metodu to uopće nije važno. Formalno, Galerkinova metoda se može primijeniti uvijek, pa i u slučaju nelinearnih rubnih problema. Tada sustav jednadžbi koji dobijemo nije više linearan.