Numericke metode u fizici, Skripte' predlog Matematički metod numeričke analize
gordana-zivkovic
gordana-zivkovic

Numericke metode u fizici, Skripte' predlog Matematički metod numeričke analize

70 str.
8broj poseta
Opis
Matematicki metod numericke analize
20 poeni
poeni preuzimanja potrebni da se preuzme
ovaj dokument
preuzmi dokument
pregled3 str. / 70
ovo je samo pregled
3 prikazano na 70 str.
ovo je samo pregled
3 prikazano na 70 str.
ovo je samo pregled
3 prikazano na 70 str.
ovo je samo pregled
3 prikazano na 70 str.
Microsoft Word - Predavanja3.doc

1

Uvodne napomene

- šta je fizika i koja je njena uloga

a) na osnovu datog stanja fizičkog sistema rekonstruisati stanje u nekom

prethodnom trenutku

b) Na osnovu datog stanja fizičkog sistema izvršiti predikciju u kakvom će

stanju sistem biti posle nekog vremena

- tipične primene komjutera u fizici:

a) Administrativno-kancelarijski poslovi

(pisanje tekstova, crtanje grafika, upoređivanje podataka, pretraživanje po

raznim ključevima...)

b) Akvizicija podataka i kontrola eksperimenta

(praćenje eksperimenta i povratnom spregom održavamo sistem u zadatom

stanju)

c) Baze podataka i obrada podataka

(računanje srednjih vrednosti, interpolacija, ekstrapolacija, fitovanje

parametara, traženje kritičnih parametra, faznih prelaza, rešavanje jednačina

koje se ne mogu analitički rešiti... Internet kao baza podataka.)

d) Kompjuterske simulacije!

- opšte napomene o važnosti numeričkih metoda, brzini izvršavanja algoritama

i programskim jezicima

- fortran - velika baza programa koji se koriste u nauci

- c - prenosiv - Unix i NT - prirodni jezici

- Pascal - brzo razvijanje programa

- Visual C, Visual Basic, Delphi

- kurs: Numerical Recipes – second edition

- Web – rudjer – Katedra za nuklearnu fiziku – Numeričke metode

2

1. Reprezentacija brojeva, opseg i tačnost

- reprezentacija brojeva u kompjuteru: dva osnovna prototipa:

brojevi sa nepokretnim i pokretnim zarezom

- celobrojne varijable

tip bita opseg

byte 8 0..255

shortint 8 -128..127

word 16 0..65535

integer 16 -32768..32767

dword 32 0..4294967296

longint 32 -2.1⋅109..2.1⋅109

- brojevi u pokretnom zarezu s⋅M⋅Be-E

s - znak

M – mantisa (normalizovana 0.5 <= M < 1)

B – osnova reprezentacije (obično 2)

E – bias (pomeranje, korekcija) eksponenta

- mašinska tačnost εm – najmanji broj koji dodat na 1 proizvodi rezultat različit

od 1 (zavisi od broja bitova mantise)

- najmanji broj koji može da se reprezentuje razlikuje se od εm

(zavisi od broja bitova eksponenta)

- standardni tipovi realnih promenjljivih u pokretnom zarezu

tip byte-a opseg εm sig. cifara single 4 10

-39 ..10

38 6⋅10-8 7-8

real 6 10 -39

..10 38

9⋅10-13 11-12 double 8 10

-324 ..10

308 10

-16 15-16

extended 10 10 -4932

..10 4932

5⋅10-20 19-20

- mašinska tačnost dovodi do grešaka zaokruživanja (roundoff errors)

- posle N uzastopnih operacija, tipična greška iznosi mNε

znak mantisa (23 bita) eksponent

3

- greške metoda (truncation error) potiču od nesavršenosti algoritama

upotrebljenih za rešavanje datog problema (razlika izmedju pravog-tačnog

rešenja i rezultata izračunavanja)

- primer: numerička integracija – kompjuter sa beskonačnom tačnošću –

rezultat izračunavanja nije beskonačno tačan

- nepovoljni slučajevi

a) može da se desi da se greška akumulira samo u jednom pravcu, Nεm

b) specifičnosti problema mogu da uzrokuju velike greške zaokruživanja

c) loši algoritmi mogu da dovedu do velikih grešaka zaokruživanja

Primer:

Sabrati 10000 puta 10 -10

(rezultat 10 -6

)

single 9.999342⋅10-7 rel. greška 7⋅10-5

real 9.99........⋅10-7 rel. greška 2⋅10-9

double rel. greška 7⋅10-14

4

2. Klasifikacija grešake, tačnost i stabilnost algoritama

- klasifikacija grešaka:

a) inherentne greške (matematiški opis problema nije egzaktan, ulazni podaci

pogrešni i slično)

b) greške zaokruživanja (roundoff errors, zavise od mašinske tačnosti)

c) greške metoda (algoritam koji se koristi ne opisuje dati problem egzaktno,

ili je potreban neprihvatljivo veliki broj operacija)

- obično ove tri vrste grešaka ne interferiraju jedne sa drugima

- u nekim slučajevima atraktivni algoritmi (brži, prostiji) za rešavanje nekog

problema mogu se pokazati nestabilni (ne konvergiraju korektnom rešenju)

- primer (veštački), zlatna sredina 61803399.0 2

15 =

− =Φ

- pretpostavimo da treba da se izračuna Φn, pokazuje se da je zadovoljena

rekurentna relacija Φn+1=Φn-1-Φn, uz Φ0=1 i Φ1=Φ

- prednost – oduzimanje umesto množenja

- pokazuje se da ovaj algoritam nije stabilan, i da u single preciznosti već Φ20

ima grešku od 100%

Primer:

Rešavanje kvadratne jednačine ax 2 +bx+c=0

a

acbb x

2

42

2/1

−±− = , ako je b2<<4ac, velike greške zaokruživanja

x 2 -800x+1=0, digitron sa 12 sigurnih cifara

digitron tačno rešenje rel. greška

x1=799.99875 799.99875 0

x2=1.25⋅10 -3

1.2500019531⋅10-3 2⋅10-6

Pravi način:

( ) q

c x

a

q xacbb

b q ==−+−= 21

2 4 2

)sgn(

5

iteracija greška

single 20 100%

real 32

double 41

extended 52

- rang algoritma – pokazuje kako raste broj operacija koje treba izvršiti u

funkciji složenosti problema

3. Sistemi linearnih algebarskih jednačina (2.0)

- sistemi oblika

MNMNMMM

NN

NN

NN

bxaxaxaxa

bxaxaxaxa

bxaxaxaxa

bxaxaxaxa

=++++

=++++

=++++

=++++

L

LL

L

L

L

332211

33333232131

22323222121

11313212111

- u matričnoj formi bxA rr

= , gde je

   

   

=

   

   

=

3

2

1

21

22221

11211

b

b

b

aaa

aaa

aaa

MNMM

N

N

L

L

L

L

L

bA

- M jednačina sa N nepoznatih

- M<N – principijelno više rešenja (može i da ih ne bude) (degenerisan)

- M>N – principijelno nema rešenja (mada može da ih bude) (predeterminisan)

- M=N – postoji dobra šansa da se sistem može rešiti

- problemi : jedan ili više redova ili kolona su linearna kombinacija ostalih

redova ili kolona – singularna matrica – det=0

- ovaj efekat se zove degeneracija po redovima ili kolonama (kod kvadratnih

matrica jedno povlači drugo)

- dodatni problemi kod nesingularnih matrica

6

a) neka od jednačina može biti vrlo blizu linearne kombinacije ostalih, pa

zbog grešaka odbacivanja i zaokruživanja, algoritam može pasti.

(Slabo uslovljeni sistemi, determinanta blizu 0).

b) akumuliranje greške računanja za veliko N može pokvariti rešenje.

Numerička procedura nije pala algoritamski, već je dobijeno rešenje

dosta različito od pravog (tačnog) rešenja.

- ako je N do 50, može da se radi u single precision, preko toga u double, sve

negde do N nekoliko stotina.

- ako je N nekoliko hiljada, problem postaje procesorsko vreme

4. Gauss-Jordan-ova eliminacija (2.1)

- pomoćni iskaz: zamena bilo kojeg reda matrice A sa linearnom kombinacijom

njega i ostalih redova matice A ne utiče na rešenje, ako se ista procedura uradi

i sa redovima vektora b

- ideja: sistem bxA rr

= transformisati u sistem 'bxI rr

= , a odatle sledi da je

' bx rr

=

- algoritam:

1) prva jednačina se podeli sa 11a (odnosno 1 '

11 =a )

2) prva jednačina se pomnoži sa 1ia i oduzme od i-te jednačine (i<>1)

(posle ovoga svi elementi u prvoj koloni osim 11a su jednaki 0)

3) druga jednačina se podeli sa 22a (odnosno 1 '

22 =a )

4) (druga jednačina se pomnoži sa ia2 i oduzme od i-te jednačine (i<>2)

(posle ovoga svi elementi u drugoj koloni osim 22a su jednaki 0)

5) postupak se nastavlja do kraja

6) na kraju matrica A postaje 1, a vektor b predstavlja rešenje

- problem nastaje ako je element 0=iia

- pivotiranje – pomoćni iskazi

7

a) zamena dva reda matrice A i odgovarajućih redova vektora b ne utiče na

rešenje

b) zamena dve kolone matrice A i odgovarajućih redova vektora x ne menja

rešenje (treba voditi evidenciju o zameni, da bi mogao da se restaurira

pravilan poredak ix

- pivotiranje – zamena redova (parcijalno pivotiranje) ili redova i kolona (puno

pivotiranje), tako da se željeni element dovede na dijagonalu

- kod punog pivotiranja mora da se vodi evidencija o zameni kolona, tako da

kasnije može da se restaurira pravilan poredak elemenata vektora x

- pivot element se obično bira tako da bude najveći koeficijent ija , ali je

poželjnije da se pre toga ceo sistem skalira, tako da u jednoj jednačini najveći

koeficijent bude 1 (implicitno pivotiranje)

- Gaus-Jordan-ova eliminacija nije stabilan algoritam ako nije uključeno

pivotiranje

- traženje inverzne matrice AA-1=1

- algoritam reda N3

5. Gauss-ova eliminacija sa povratnom zamenom (2.2)

- slično kao kod Gauss-Jordan-ove eliminacije, samo što se u prvom prolazu

eliminišu jednačine koje su ispod tekuće jednačine

- posle prvog prolaza, dobijamo sistem

    

    

=

   

   

    

    

'

4

'

3

'

2

'

1

4

3

2

1

'

44

'

34

'

33

'

24

'

23

'

22

'

14

'

13

'

12

'

11

000

00

0

b

b

b

b

x

x

x

x

a

aa

aaa

aaaa

- povratna zamena – očigledno je x4

'

44

'

44 / abx =

- kada znamo x4, onda je x3

8

[ ]'344'3' 33

3

1 axb

a x −=

- odnosno, ako produžimo dalje

 

  

 −= ∑

+=

N

ij

jiji

ii

i xab a

x 1

''

'

1

- i Gauss-ov algoritam bez pivotiranja predstavlja nestabilan algoritam

- algoritam reda N3

6. LU dekompozicija (2.3)

- pretpostavimo da matricu A možemo da napišemo u obliku

LU=A

- L i U su triangularne matrice takve da je

  

  

=

  

  

  

  

333231

232221

131211

33

2322

131211

333231

2221

11

00

00

00

aaa

aaa

aaa

β

ββ

βββ

ααα

αα

α

- ovu dekompoziciju koristimo da rešimo sistem

Ax=(LU) x=L(Ux)=b

- prvo rešavamo Ly=b, a zatim Ux=y

- prednost – kod tridiagonalnih matrica, rešavanje je znatno olakšano

- prvu jednačinu lako rešavamo zamenom unapred

Niyby

b y

i

j

jiji

ii

i ,...,3,2 1 1

1

11

1 1

= 

  

 −=

=

∑ −

=

α α

α

- drugu jednačinu rešavamo zamenom unazad

9

1,...,2,1 1

1

−−= 

  

 −=

=

∑ +=

NNixyx

y x

N

ij

jiji

ii

i

NN

N N

β β

β

- ovaj algoritam je reda N2, pod uslovom da se znaju matrice L i U

- kako odrediti matrice L i U – vidi se da je

ijji a=+ ...11βα

- broj članova zavisi od toga da li je i>j ili obrnuto

)3(:

)2(:

)1(:

2211

2211

2211

ijjjijjiji

ijjjiijiji

ijijiijiji

aji

aji

aji

=+++>

=+++=

=+++<

βαβαβα

βαβαβα

βαβαβα

L

L

L

- ukupno imamo N2 jednačina sa N2 + N nepoznatih

- stavljamo da je 1=iiα

- Crout-ov algoritam

1) postavljamo 1=iiα

2) za svako j=1,2,..,N uraditi sledeće (formule 1 i 2)

a) za svako i=1,2,..,j naći ijβ preko formule (1 i 2)

∑ −

=

−= 1

1

i

k

kjikijij a βαβ

b) za svako i=j+1,..,N (formula 3)

 

  

 −= ∑

=

1

1

1 j

k

kjikij

jj

ij a βαβ α

- pažljivim praćenjem ovog algoritma, primećujemo da su α i βna desnoj strani

uvek određeni u trenutku kada zatrebaju za račun. Sa druge strane, elementi

ija uvek se javljaju samo jednom, pa posle njihovog korišćenja, njihove

memorijske lokacije se mogu koristiti za smeštanje koeficijenata α i β.

10

   

   

44434241

34333231

24232221

14131211

βααα

ββαα

βββα

ββββ

- pivotiranje je apsolutno neophodno za stabilnost Crout-ovog algoritma, ali je

samo parcijalno pivotiranje efikasno upotrebljivo. Medjutim, i to je dovoljno

da obezbedi stabilnost.

- algoritam red N3

- determinanta matrice ∏ =

= N

j

jjb 1

det

7. Iterativno poboljšanje rešenja sistema linearnih jednačina (2.5)

- jasno je da se ne može dobiti rešenje sistema linearnih jednačina tačnije od

mašinske tačnosti upotrebljenog kompjutera

- za veliko N, teško je dobiti rešenje čija je tačnost uporediva sa mašinskom

tačnošću

- iterativno poboljšanje

- neka je x tačno rešenje

Ax=b

- naše rešenje – približno x+δx, gde je δx nepoznato

A(x+δx)=b+δb

- odnosno

A⋅δx=δb

- ako imamo LU dekompoziciju, δx se lako nalazi

- red algoritma N2, pa se preporučuje (već je potrošeno N3 operacija)

11

- pretpostavka je da je LU dekompozicija apsolutno tačna, što ne mora biti

slučaj

- neka je B0 aproksimativno A -1

, tako da je AB0 približno 1

- definišimo rezidualnu matricu R=1-B0A, koja je mala

- tada je B0A=1-R

- napravimo sledeću formalnu manipulaciju

0

2

0

1

0

1

00

1

0

1

0

1

0

11 B...)RR(1BR)(1BA)(BB)B(A)B(BAA ⋅+++=⋅−=⋅⋅=⋅⋅=⋅⋅= −−−−−−−

- označimo n-tu parcijalnu sumu sa Bn

0

n2

n B)R...RR(1B ⋅++++=

- pa je xn=Bn⋅b

Primer:

  

  

=

  

  

1

1

1

300

020

001

x , očigledno je da je   

  

=

3/1

2/1

1

x ,   

  

=−

3/100

02/10

001 1

A

pođimo od ‘približnog’ rešenja   

  

=+=

1

1

1 ' xxx δ

  

  

=

  

  

  

  

=

  

  

  

  

  

  

=−+=

2

1

0

1

1

1

3

2

1

1

1

1

1

1

1

300

020

001

)( bxxAb δδ

  

  

=

  

  

  

  

== −

3/2

2/1

1

2

1

0

3/100

02/10

001 1 bAx δδ , odnosno

  

  

=

  

  

  

  

=

3/1

2/1

1

3/2

2/1

1

1

1

1

x

12

8. Sherman-Morrison-ova formula (2.7)

- pretpostavimo da imamo smo dobili inverznu matricu A-1 neke matrice A

- pretpostavimo da se matrica B malo razlikuje od matrice A (jedan element,

jedan red, jedna kolona)

- može li se dobiti inverzna matrica matrice B iz inverzne matrice A-1

vuAB ⊗+=

- gde je

[ ]N

N

vvv

u

u

u

L M

21

2

1

   

   

=⊗ vu

- ima smisla ako je ili u ili v vektor koji ima samo jednu komponentu različitu

od 0 (u tom slučaju, matrice B i A se razlikuju u jednom redu ili jednoj

koloni)

λ

λλ

+

⋅⊗⋅ −=

−+−⋅⊗⋅−=

⋅+⊗⋅⋅⊗⋅+⊗⋅−=

⋅⊗⋅+=⊗⋅+⋅=⊗+

−− −

−−−

−−−−

−−−−−−

1

)()(

...)1(

...)(

)())(()(

11 1

2111

1111

11

AvuA A

AvuAA

AvuAvuAvuA1

AvuA1vuA1AvuA 1111

- gde je

uAv ⋅⋅= −1λ

- operativno

zvvAwuAz 1 ⋅=⋅=⋅= −− λT)( 1

- pa je

λ+ ⊗

−= −− 1

1 wzAB 1

13

Primer:

[ ]100 0

1

1

100

010

001

000

100

100

100

010

001

100

110

101

  

  

+

  

  

=

  

  

+

  

  

=

  

  

pa je

  

  

  

  

==⋅==⋅=

  

  

=

  

  

= −−

1

0

0

0

1

1

)(

1

0

0

0

1

1 1 λvvAwuuAzvu 1 T

  

  

= +

⊗ −=

  

  

=⊗=⊗ −−

100

110

101

1 000

100

100 1

λ wz

ABvuwz 1

Primer:

[ ] [ ]4320

0

0

0

1

0000

0000

0000

4320

4000

0300

0020

0001

4000

0300

0020

4321

   

   

+=

   

   

+

   

   

=

   

   

   

   

=

   

   

=

   

   

=−

4

3

2

0

0

0

0

1

4/1000

03/100

002/10

0001

1 vuA

0

1

1

1

0

)( 11 =⋅=⋅=

   

   

=⋅==⋅= −− uvzvvAwuuAz λT

   

   

 −−−

=

   

   

=⊗ −

4/1000

03/100

002/10

1111

0000

0000

0000

1110

1Bwz

14

9. Interpolacija i ekstrapolacija funkcija (3.0)

- pretpostavimo da znamo vrednosti neke f-je u N tačaka, a ne znamo

analitički oblik f-je

- (x1,y1), (x2,y2),...., (xN,yN) – sortiramo tako da je x1≤x2≤...≤xN

- primer: imamo rezultate N merenja vrednosti neke veličine y, za N

vrednosti nekog parametra x

- problem – naći vrednosti f-je y za proizvoljno x

- ako ne x1 ≤ x ≤ xN, onda je u pitanju interpolacija, u suprotnom

ekstrapolacija

- sa ekstrapolacijom treba biti obazriv, može da da potpuno pogrešne

rezultate (koji značajno odstupaju od merene veličine y)

- ako znamo nešto o f-ji y, onda možemo da primenimo fitovanje (mnogo

bolje, pošto merene veličine obično imaju nezanemarljive greške merenja)

- za dati set podataka, interpolaciona f-ja prolazi kroz tačke (xi,yi), a

fitovana ne mora

- obično ne znamo mnogo o analitičkom obliku f-je y (ili je ona mnogo

komplikovana), pa datu (nepoznatu) f-ju aproksimiramo drugom f-jom

- obično se koristi aproksimacija polinomima, racionalnim f-jama i spline

funkcijama

- pokazuje se da je zbog grešaka računanja bolje računati vrednosti f-je za

svaku traženu vrednost x, nego naći parametre f-je, pa je onda računati za

svako traženo x

10. Lagrange-ova interpolaciona formula i Neville-ov algoritam (3.0)

- kroz dve tačke možemo da provučemo liniju (kriva prvog reda), kroz tri

tačke kvadratnu f-ju (2 reda)

- N tačaka, kriva N-tog reda

15

- ako imamo N parova vrednosti (x1,y1), (x2,y2),...., (xN,yN), onda kroz te

tačke možemo provući polinom N-1 reda y=c0+c1x+...+cN-1x N-1

- Lagrange-ova formula

N

NNNN

N

N

N

N

N

y xxxxxx

xxxxxx

y xxxxxx

xxxxxx y

xxxxxx

xxxxxx xP

))...()((

))...()(( ...

))...()((

))...()((

))...()((

))...()(( )(

121

121

2

23212

31

1

13121

32

−−−

−−− ++

+ −−−

−−− +

−−−

−−− =

- u tačkama u kojima je x=x1,x2,...,xN ostaje samo jedan član <>0 (=yi)

- može direktno da se koristi, ali je bolje primeniti Neville-ov algoritam

- neka je Pi polinom nultog reda (konstanta) koji prolazi kroz tačke (xi,yi)

- Pi=yi

- neka je Pi,i+1 polinom prvog reda koji prolazi kroz tačke (xi,yi) i (xi+1,yi+1)

)(

)()(

)(

)(

)(

)(

)(

)(

)(

)(

1

11

1

11

1 1

11

1 1,

+

++

+ ++

+ +

++

+ +

−+− =

= −

− +

− =

− +

− =

ii

iiii

i

ii

i i

ii

i i

ii

i i

ii

i ii

xx

PxxPxx

P xx

xx P

xx

xx y

xx

xx y

xx

xx P

- lako se pokazuje da važi

)(

)()( ))...(2)(1()1)...(1( ))...(1(

mii

miiiimiiimi

miii xx

PxxPxx P

+

+++−+++ ++ −

−+− =

- tok izračunavanja

444

34

234333

123423

123222

12

111

:

:

Pyx

P

PPyx

PP

PPyx

P

Pyx

=

=

=

=

- ako su nam potrebni koeficijenti interpolacionog polinoma

y=c0+c1x+...+cNx N

- gde ci zadovoljavaju

16

   

   

=

   

   

    

    

NN N

NN

N

N

y

y

y

c

c

c

xx

xx

xx

LL

L

L

L

L

1

0

1

0

11

00

1

1

1

- treba rešiti ovaj sistem jednačina

- Vandermond-ova matrica

- slabo uslovljena (za veliko N velike greške zaokruživanja, koristiti double

precision)

- primer: za dvadeset tačaka x od 1 do 20, javlja se 2020≈1026

- posebni algoritmi

11. Interpolacija kubnim spline-om (3.3)

- neka je data tabelirana f-ja yi=y(xi), i=1,2,...,N

- posmatrajmo samo interval između xi i xi+1

- ako interpoliramo linearnom f-jom samo ovaj interval, onda je

1++= jj ByAyy

- gde je

jj

j

jj

j

xx

xx AB

xx

xx A

− =−=

− =

++

+

11

1 1 (1)

- ove jednačine su specijalni slučaj Langrange-ove jednačine

- pošto je gornja jednačina linearna po x, to su njen prvi i drugi izvod

0''' 1

1 = −

− =

+

+ y

xx

yy y

jj

jj

- drugi izvod je nula svuda osim u tačkama xj, gde je nedefinisan

- interpolacijom kubnim spline-om dobijamo f-ju kojoj je prvi izvod gladak

(nema prelomnih tačaka) a drugi izvod neprekidan

- pretpostavimo da pored f-je yi znamo i vrednosti drugog izvoda u y’’ (što

nije tačno)

17

- želimo da se drugi izvod y’’ u intervalu izmedju xj i xj+1 menja linearno

počev od vrednosti ''

jy do vrednosti ''

1+jy

''

1

''''

++= jj ByAyy

- lako se pokazuje da onda interpolacioni polinom mora da ima vrednost

''

1

''

1 ++ +++= jjjj DyCyByAyy

- gde su A i B vrednosti jednačinama (1), a C i D

2

1

32

1

3 ))(( 6

1 ))((

6

1 jjjj xxBBDxxAAC −−=−−= ++

- zavisnost trećeg reda od x

- ako se interpolacioni polinom dva puta diferencira, dobija se

''

1

''

2

2

''

11

2 ''

1

2

1

1 )(

6

13 )(

6

13

+

+++ +

+

+=

− −

+− −

− −

− =

jj

jjjjjj

jj

jj

ByAy dx

yd

yxx B

yxx A

xx

yy

dx

dy

(2)

- drugi izvod je baš ono što smo tražili

- jedini problem je što ne znamo y’’

- ideja je da se traži neprekidnost prvog izvoda u tačkama xj, odnosno da

vrednosti prvog izvoda u tački xj računata iz intervala (xj-1,xj) i, posebno, iz

intervala (xj,xj+1), budu jednake

- tada iz (2) dobijamo

1

1

1

1''

1

1''11''

1

1

636 −

+

+ +

+−+ −

− −

− =

− +

− +

jj

jj

jj

jj

j

jj

j

jj

j

jj

xx

yy

xx

yy y

xx y

xx y

xx

- N-2 linearne jednačine za N nepoznatih y’’

- dva načina prevazilaženja ovog problema

a) u graničnim tačkama se pretpostavi da je 0 ''

1 =y i 0 '' =Ny

(prirodni kubni spline)

b) postavi se da ''

1y i ''

Ny imaju neke vrednosti izračunate preko (2),

tako da prvi izvodi u graničnim tačkama imaju zadate vrednosti

18

- pošto je matrica tridijagonalna, rešavanje ovog sistema zahteva N

operacija

12. Elementarni algoritmi za integraciju funkcija (4.0)

- numerička integracija (kvadratura) spada u jedan od najstarijih problema

numeričke analize

- problem je izračunati integral

∫= b

a

dxxfI )(

- što je ekvivalentno rešiti diferencijalnu jednačinu

0)(),( == ayxf dx

dy

- problem: integral ne može da se izračuna u analitičkom obliku, ili je f-ja

data u tabelisanom obliku (rezultat merenje) (integral ispod pika)

- klasične formule za ekvidistantne tačke skoro neupotrebljive, mada vrlo

lepe i elegantne (muzejski primerci)

- notacija: abscisa x1,...,xN+1 ; konstantni korak h

1,...,1,00 +=+= Niihxxi

- dva tipa formula: zatvorene i otvorene, u zavisnosti da li krajnje tačke

ulaze u formulu (na krajevima f-ja ima singularitet)

19

- formule Newton-Cotesovog tipa – aproksimacija f-je polinomom

odgovarajućeg stepena, a zatim integracija tog polinoma

- Trapezno pravilo (polinom prvog reda, tačna za polinom prvog reda)

)''( 2

1

2

1 )( 321

2

1

fhOffhdxxf

x

x

+

 

 +=∫

- greška zavisi od h3 i drugog izvoda f-je f (koji ne znamo) negde izmedju

x1i x2

- Simsonovo pravilo (polinom drugog reda, tačna do polinoma trećeg reda)

Primer: Izvesti trapezno pravilo

2

12

1 1

21

2 2211 )(),(),( y

xx

xx y

xx

xx xfyxyx

− +

− =

[ ]

[ ]

)( 2

1 ))((

2

1

2

1

22 2

1

)())(( 2

1

)())(( 2

1

)()( 2

1

2

)(

121212

11211222

211211211222

21121212

21121221

12

21

21122

1

2

2

21

21

21

2112

2

21

21

21

21

21

12

21

2

21

1

2

12

1 1

21

2

2

1

2

1

2

1

2

1

yyhxxyy

yxyxyxyx

yxyxyxyxyxyx

yxyxxxyy

yxyxxxyy

xx xx

yxyx xx

xx

yy

xx

yxyxx

xx

yy

dx xx

yx

xx

yx x

xx

y

xx

y

dxy xx

xx y

xx

xx dxxf

x

x

x

x

x

x

x

x

+=−+=

−−+=

−+−+−=

−++−=

−++−−=

− −

− −−

− =

− −⋅

− =

 

  

  

  

− −

− −

  

− −

− =

 

  

− +

− =

∫ ∫

20

)( 3

1

3

4

3

1 )( )4(5321

3

1

fhOfffhdxxf

x

x

+

 

 ++=∫

- Simpsonovo 3/8 pravilo (pol. 3 reda, tačna do polinoma 3 reda)

)( 8

3

8

9

8

9

8

3 )( )4(54321

4

1

fhOffffhdxxf

x

x

+

 

 +++=∫

- Bode-ovo pravilo (pol. 4 reda, tačna za pol. do 5 reda)

)( 45

14

45

64

45

24

45

64

45

14 )( )6(754321

5

1

fhOfffffhdxxf

x

x

+

 

 ++++=∫

13. Prošireno trapezno pravilo (4.1)

- izdelimo abscisu na N podintervala, onda je po trapeznom pravilu

∫ −+++= Nx

x

NIIIdxxf

1

121 ...)(

- odnosno

)''( 2

1

2

1 ...

2

1

2

1 )( 3121

1

fhNOffhffhdxxf NN

x

x

N

+

 

 +++

 

 += −∫

- prošireno trapezno pravilo

 

  

 − +

 

 ++++=∫ − 2 3

2121

'')(

2

1 ...

2

1 )(

1 N

fab Offffhdxxf

Nx

x

N

- ova formula se koristi u mnogim algoritmima algoritmima za num. integraciju

- proširena Simsonova formula

 

  

+

 

 +++++++=∫ −− 4124321 1

3

1

3

4

3

2 ...

3

4

3

2

3

4

3

1 )(

1 N

Offfffffhdxxf Nx

x

NNN

- problem : kako odrediti broj tačaka za datu grešku

- način realizacije integracije pomoću proširenog trapeznog pravila (slika)

- S0 – integral preko dve tačke, podelimo interval na pola – S1, opet na pola – S2

21

- Si – ukupan broj tačaka 2 i +1, korak iii

hab h

22

0= −

=

) 2

()( 2

1 )(

2

1 )(

2

1 )

2 ()(

2

1

)( 2

1 )(

2

1

2

1

2

1

1111

02100

ab afhbfafhbf

ab afafhS

bfafhffhS

− ++

 

 +=

 

 + −

++=



 

 +=

 

 +=

- imajući u vidu da je h0=2h1

) 2

( 22

0

1

ab af

abS S

− +

− +=

- analogno



 

 − ++

− +

− +=



 

 − ++

− ++

 

 +

− ++=



 

 + −

++ −

++ −

++=

) 4

3() 4

( 42

) 4

3() 4

()() 2

()(

)() 4

3() 4

2() 4

()(

1

22

22

ab af

ab af

abS

ab af

ab afhbf

ab afafh

bf ab

af ab

af ab

afafhS

- u opštem slučaju

∑ −

= −

− −+

− +=

12

1

1 ) 2

)12(( 22

1 i

j iiii

ab jaf

ab SS

- metodom sukcesivnih aproksimacija prilazimo pravoj vrednosti integrala

MSSS →→→ ...10

- stati kada je ε<− −1MM SS , gde je ε zadata greška

Primer: naći 3

2

3

1

1

31

1

2 == −−

x

dxx

iteracija tačaka korak I

1 2 2 2

2 3 1 1

3 5 0.5 0.75

4 9 0.25 0.6875

5 17 0.125 0.671875

6 33 0.0625 0.667969

22

14. Rombergova formula (4.3)

- greška prilikom računanja proširenom trapeznom formulom srazmerna h 2

- neka je I prava vrednost integrala

2

22

2

11

ahSI

ahSI

=−

=−

- pa je

2

2

2

1

2

1

SahI

SahI

=−

=−

2

1

1

h

SI a

− =

 

 

 

  

 −

 

  

 −

= 

  

 −=

 

 

 

  

 −=

− −

2

1

2

2

1

2 122

1

2 12

2

1

2 2

2

22

1

1

1

1

h

h

h

h SS

I h

h SS

h

h ISh

h

SI I

- ako je h2/h1=1/2, što je slučaj kod sukcesivne primene proš. trapezne formule

12

12

3

1

3

4

4

1 1

4

1

SS

SS

I −= −

− =

- Rombergova formula

2

2

1

2

1

1

0

2

0

1

0

0

S

SS

SSS

↓↓

→→

- gde je

0

1

01

3

1

3

4 −−= iii SSS

23

- već posle prve popravke poklapa se sa proširenom Simpsonovom formulom –

tačna za polinome 2 reda

- greška od 1

1S i 1

2S srazmerna sa h 4

- ponovimo gornju proceduru i dobijamo opštu formulu

[ ]11111 12 1 −

− −

+ − −

− −= ji

j

ij

j

i

j

i SSSS

- metoda se zove Richardson-ova ekstrapolacija prema granici

- gornji primer za trapeznu formulu, primenom Rombergove formule daje tačno

rešenje već u drugoj iteraciji

15. Izračunavanje funkcija pomoću stepenih redova (5.0)

- funkcija f u okolini tačke x0 može da se predstavi stepenim redom

∑ ∞

=

−= 0

0 )()( k

k

k xxaxf

- zahtev – red bi trebao da konvergira i to po mogućstvu brzo (viši članovi brzo

teže prema 0)

- Tejlorov razvoj

∑ ∞

=

−= 0

0 0

)(

)( !

)( )(

k

k k

xx k

xf xf

Primer: naći 2)cos()sin( 00

=−=∫ ππ

xdxx

Iteracija Trapez ε Romberg ε

1 0 1 0 1

2 1.5708 2⋅10-1 2.0943 5⋅10-2

3 1.8961 5⋅10-2 1.9917 4⋅10-3

4 1.9742 10 -2

2.0002 9⋅10-5

5 1.9936 3⋅10-3 1.9999 2⋅10-6

6 1.9984 8⋅10-4 2.0000 10-8

24

- ako je x0=0 – Tejlor-Maklorenov razvoj

- oblast kovergencije definiše vrednosti x za koje red konvergira

- neke standardne f-je

11 !

)1)...(1( 1...

!2

)1( 1)1(

11 1

)1(... 32

)1ln(

)!2( )1(...

!4!2 1)cos(

)!12( )1(...

!5!3 )sin(

! ...

!3!2 1

0

2

0

132

0

242

0

1253

0

32

<<− +−−

+=+ −

++=+

<<− +

−=−+−=+

∞<<−∞−=−+−=

∞<<−∞ +

−=−+−=

∞<<−∞=++++=

=

=

+

=

=

+

=

xx k

k xxx

x k

xxx xx

x k

xxx x

x k

xxx xx

x k

xxx xe

k

k

k

k k

k

k k

k

k k

k

k x

ααααα αα

- obratiti pažnju na način izračunavanja članova: 3!2!3

23 xxx ⋅=

- pre izračunavanja pokušati da se x smanji pogodnim operacijama (ako može)

- sin i cos svesti na opseg (0,Pi)

- eksponencijalna f-ja i logaritam – iskoristiti ln(ab)=ln(a)+ln(b)

- ubrzavanje konvergencije stepenih redova

- Aitken-ov δ2 proces pretpostavka da red konvergira – znamo tri sukcesivne

parcijalne sume – ekstrapoliramo za ∞→n

( ) 11

2

1

1

'

2 −+

+ + +−

− −=

nnn

nn

nn SSS

SS SS

- za alternativne redove može da se primeni Euler-ova transformacija

- primena

a) izračunavanje f-ja

b) numeričko integraljenje

∫ ∑∑∫ ∞

=

+ ∞

=

−=−= 0

1

0

0

0 )()()( k

kk

k

k

k xx k

a dxxxadxxf

c) numeričko diferenciranje

25

∑∑ ∞

=

− ∞

=

−=−= 0

1

0

0

0 )()()( k

k

k

k

k

k xxkaxxa dx

d xf

dx

d

16. Izračunavanje funkcija pomoću Chebyshev-ljevih polinoma

(5.8)

- Chebyshev-ljevi polinomi reda n definisani su kao

))arccos(cos()( xnxTn =

- liči na čistu trigoometrijsku f-ju, ali se pokazuje da je

)()(2)(

34)(

12)(

)(

1)(

11

3

3

2

2

1

0

xTxxTxT

xxxT

xxT

xxT

xT

nnn −+ −=

−=

−=

=

=

L

nema postavljenih komentara
ovo je samo pregled
3 prikazano na 70 str.