Skozi grede > Ropotarnica

Numerični algoritmi

Interpolacija funkcijReševanje enačbSistem linearnih enačbMatrike in vektorjiOdvod in integralEnačba rastiGibalna enačbaAdvekcijska enačbaDifuzijska enačbaValovna enačbaPotencialna enačbaKompleksna difuzijska enačbaKompleksna potencialna enačbaSpektralna analiza – Slučajna števila

Računalnik je neprekosljivo orodje za pisanje, risanje in računanje (ter še za mnogo drugega). Zlasti računanje je vsakdanji kruh matematikom, fizikom in tehnikom. Poglejmo, kako rešujemo tipične številske probleme v fiziki! Omejili se bomo na preproste algoritme. Marsikod obstajajo boljši pristopi, celo mnogo boljši, a so temu primerno tudi bolj zamotani.

Interpolacija funkcij

Funkcija u(x), ki je opisana s tabelo, je poznana zgolj v posameznih točkah. Zgled so logaritemske ali trigonometrične tablice. V vmesnih točkah računamo funkcijske vrednosti z interpolacijo.

Linearna interpolacija.

Med dvema točkama u1 = u(x1) in u2 = u(x1 + h) funkcijo aproksimiramo s premico. V točki x1 + θh, 0 ≤ θ ≤ 1, ima funkcija potem vrednost u:

u = u1 + θ(u2u1) .

Parabolična interpolacija.

Bolj natančna je aproksimacija s parabolo skozi tri zaporedne ekvidistantne točke:

u = u1 + θ(u2u1) + (1/2)θ(θ − 1)(u3 − 2u2 + u1) .

Reševanje enačb

Kadar kakšne enačbe f(x) = 0 ne znamo ali ne zmoremo rešiti algebraično, s simboli, jo rešujemo numerično, s števili.

Iteracija.

Če lahko funkcijo izrazimo kot f(x) = x + R(x), dobimo:

x = −R(x) .

V desno stran vstavimo primeren približek x0 in izračunamo levo stran – novi približek x1. Tako nadaljujemo. Metoda konvergira, če je v okolici rešitve |R'(x)| ≤ 1.

Tangentna metoda.

Zaporedne približke k rešitvi določamo tudi s pomočjo tangent na graf funkcije. Tangenta v točki f(x0) seka koordinatno os v točki x1 in to je tudi novi približek:

x1 = x0

f(x0)

f'(x0)

.

Sekantna metoda.

Iteracija in tangentna metoda ne delujeta v vseh primerih. Zmeraj pa deluje naslednja metoda. Z grobim tabeliranjem najdemo dve vrednosti a in b, pri katerih ima funkcija nasprotna predznaka: rešitev funkcije leži tedaj nekje na intervalu [a, b]. Premica skozi navedeni točki seka koordinatno os v točki c. Odvisno od predznaka funkcije v tej točki je novi interval [a, c] ali [c, b]:

c = af(a)

ba

f(b) − f(a)

.

Sistem linearnih enačb

Sistem n × n linearnih enačb je v matrični obliki zapisan kot A · x = b. Njegova formalna rešitev je x = A−1 · b. Numerično jo določimo takole.

Eliminacija.

Najprej k matriki A na desni strani prilepimo stolpec b in tako dobimo razširjeno matriko A|b koeficientov. Potem:

  1. Najdemo "delovno" vrstico z (absolutno) največjim vodilnim elementom in jo postavimo na vrh.
  2. Delovno vrstico delimo z njenim prvim elementom.
  3. Od vsake naslednje vrstice odštejemo delovno vrstico, pomnoženo s prvim elementom te vrstice.
  4. Pokrijemo prvo vrstico in prvi stolpec ter nadaljujemo s preostankom po istem postopku, le da odštevamo od vsake naslednje in od vsake predhodne vrstice.
  5. Ponovimo postopek od spodaj navzgor.

Dobimo enotno razširjeno matriko I|d, to je sistem I · x = d, ki je že iskana rešitev.

Iteracija.

Pri eliminaciji se zaokrožitvene napake akumulirajo. Kadar je matrika velika, postane rešitev neuporabna. Takrat pomaga iterativna metoda – taka, kot pri iskanju korena funkcije. Matriko zapišemo v obliki A = D + R, pri čemer vsebuje prva matrika diagonalo, druga pa preostanek. Tako dobimo enačbo D · x = bR · x. V desno stran vstavimo primeren približek x0 in izračunamo levi vektor – z diagonalnimi koeficienti pomnoženi novi približek x1. Eksplicitno velja:

xi1 =

1

aii

(bi

 

ji

aijxj0) .

Metoda konvergira, če je vsak diagonalni element matrike po absolutni vrednosti večji od vsote absolutnih vrednosti preostalih elementov v vrstici. Po želji v desno stran sproti vstavljamo že izračunane nove komponente. Lahko pa celo dobljeni novi približek xi1 poimenujemo kot "provizoričnega" in iz njega izračunamo "pravega",

xi2 = xi0 + ω(xi1xi0), 1 ≤ ω < 2 .

Matrike in vektorji

Inverzna matrika.

Sistem linearnih enačb lahko rešimo tudi z izračunom inverzne matrike. In obratno: računanje inverzne matrike prevedemo na reševanje sistema linearnih enačb. Za inverzno matriko X k matriki A velja, po definiciji, A · X = I. Za vsak (neznani) stolpec x inverzne matrike in ustrezni (znani) stolpec i enotne matrike torej velja A · x = i, kar rešujemo z eliminacijo ali iteracijo. Pri eliminaciji lahko določimo celo vse stolpce hkrati: tvorimo razširjeno matriko A|I in jo z eliminacijo preoblikujemo v I|X.

Lastne vrednosti in vektorji.

Vsako simetrično matriko Aij = Aji lahko preoblikujemo v diagonalno matriko in s tem najdemo njene lastne vektorje in lastne vrednosti. Matriko je treba "le" obdelati s primernimi rotacijskimi matrikami.

Rotacijska matrika je enotna matrika, ki ima diagonalna elementa Rpp = Rqq = cos φ = c ter izvendiagonalna elementa Rpq = −Rqp = sin φ = s. Označimo jo kot Rpq. Transformacija Rpq · A · RTpq izdela matriko A', ki je enaka izvorni matriki s spremenjenima vrsticama p in q ter stolpcema p in q. Izbrati želimo takšno rotacijsko matriko, torej takšen zasuk φ oziroma takšni vrednosti c in s, da bo element Apq postavljen na nič. Iskani zasuk je:

cot 2φ

AqqApp

2Apq

.

S tem so določene količine

c = cos φ
s = sin φ
t = s/c
τ = s/(1+c) .

Preoblikovana matrika ima naslednje spremenjene elemente (pri čemer rp, rq):

A'pq = 0
A'pp = ApptApq
A'qq = AqqtApq
A'rp = Arps(Arq + τArp)
A'rq = Arps(ArpτArq) .

Diagonalizacija poteka takole. V izvorni matriki A poiščemo največji element Apq nad diagonalo in izračunamo novo matriko A', ki ima ustrezen element postavljen na nič. Pri tem se nekateri preostali elementi spremenijo. Postopek ponavljamo na novi matriki, dokler ta ne postane diagonalna. Tako dobimo lastne vrednosti; to so diagonalni elementi. Lastne vektorje pa potem določimo iz definicijske enačbe A · x = λx, ki jo zapišemo v obliki (A − λI) · x = 0. Sistem rešimo za vsak λ na že znani način.

Odvod in integral

Enostranska in centralna razlika.

Odvod funkcije aproksimiramo z ustrezno diferenco, enostransko ali centralno:

du

dx

u(x + h) − u(x)

h

du

dx

u(x + h) − u(xh)

2h

d2u

dx2

u(x + h) − 2u(x) + u(xh)

h2

.

Trapezna in parabolična formula.

Ploščino pod krivuljo računamo tako, da krivuljo na vsakem primerno ozkem intervalu nadomestimo s premico ali parabolo, nakar ploščinske prispevke z vseh intervalov seštejemo:

x+h

x

u(x) dx

u(x) + u(x+h)

2

h

x+2h

x

u(x) dx

u(x) + 4u(x+h) + u(x+2h)

3

h .

Enačba rasti

Enačbo rasti (za spreminjanje mase vode v rezervoarju z dotoki in odtoki)

dm

dt

= f(t, m)

Tetivna metoda.

najenostavneje računamo po tetivni metodi, ki pa je le malo natančna:

m(t + dt) = m + f(t, m)dt .

Dvotetivna metoda.

Boljša je dvotetivna metoda:

m(t + dt) = m + (

k1

2

+

k2

2

)dt
k1 = f(t, m)
k2 = f(t + dt, m + k1dt) .

Gibalna enačba

Gibalna enačba (za gibanje točkastega telesa pod vplivom sile)

d2s

dt2

= f(t, s,

ds

dt

)

je ekvivalentna sistemu dveh sklopljenih enačb

ds

dt

= v

dv

dt

= f(t, s, v)

Prepletna metoda.

in tako jo tudi rešujemo. Ob času t je začetna lega s0 in ob času t + dt / 2 "začetna" hitrost v1/2. Nove vrednosti dt kasneje so:

s1 = s0 + v1/2dt
v1+1/2 = v1/2 + f1dt .

Če "začetne" hitrosti ne poznamo, jo določimo iz prave začetne hitrosti v0 s posebnim korakom: v1/2 = v0 + f0dt / 2.

Advekcijska enačba

Advekcijsko enačbo (za širjenje koncentracije primesi s snovnim tokom)

Q

t

= −c

Q

x

Shema FTUS.

aproksimiramo v točkah i z diferencami naprej v času in proti toku v prostoru. Shema ohranja pozitivnost in integral advektirane količine in je stabilna, če je časovni korak dovolj kratek:

Q1i = Qir(Qi+1Qi)
r = cdt / dx ≤ 1 .

Implicitna shema.

Če prostorski gradient aproksimiramo s centralno razliko, enkrat ob času t in enkrat ob času t + dt, ter uporabimo njuno povprečje, dobimo implicitno shemo

rQi−11 + Qi1 + rQi+11 = Qi−1 + QiQi+1
r = cdt / 2dx .

To je sistem enačb s tridiagonalno matriko aixi−1 + bixi + cixi+1 = di, a1 = 0, cn = 0. Rešimo ga z eliminacijo, ki je prilagojena za tovrstno obliko matrike. Najprej za člene od 2 do n po vrsti naprej izračunamo

m = ai / bi−1
bi = bimci−1
di = dimdi−1 ,

nato pa za člene od n−1 do 1 po vrsti nazaj

xi = (dicixi+1) / bi .

Shema je stabilna za poljuben časovni korak.

Difuzijska enačba

Difuzijsko enačbo (za širjenje koncentracije primesi v mirujoči snovi)

Q

t

= D

2Q

x2

Shema FTCS.

rešujemo v točkah i z diferencami naprej v času in centralno v prostoru ter s primerno kratkim časovnim korakom, da zagotovimo stabilnost:

Q1i = Qi + r(Qi+1 − 2 Qi + Qi−1)
r = Ddt / dx2 ≤ 1 / 2.

Implicitna shema.

Po zgledu advekcijske enačbe tvorimo implicitno shemo kot

rQi−11 + (1 + 2r)Qi1rQi+11 = rQi−1 + (1 − 2r)Qi + rQi+1
r = Ddt / 2dx2 .

Shema je stabilna za vsak časovni korak.

Valovna enačba

Valovna enačba (za širjenje snovnih ali elektromagnetnih valov)

2E

t2

=

1

c2

2E

x2

se zapiše kot sistem dveh enačb

B

t

= −

E

x

E

t

= −c2

B

x

.

Prepletna shema.

Naj bosta začetni polji v dveh naborih točk E0i in B0i+1/2; nabora točk sta medsebojno zamaknjena za interval dx / 2. V časovnem intervalu dt se električno polje spremeni v

E1i = E0ic2 (dt / dx) (B0i+1/2B0i−1/2) ,

nakar se spremeni še magnetno polje v

B1i+1/2 = B0i+1/2 − (dt / dx) (E1i+1E1i) .

Na robovih je treba upoštevati primerne robne pogoje. Da bo rešitev stabilna, moramo uporabljati dovolj kratke časovne korake: cdt / dx ≤ 1.

Potencialna enačba

Potencialno enačbo (za deformacijsko polje snovi ali za elektrostatično polje v kondenzatorju)

2ϕ

x2

+

2ϕ

y2

= 0

Relaksacija.

aproksimiramo v točkah i, j s centralnimi diferencami v prostoru ter rešujemo z relaksacijo:

ϕij1 =

1

4

(ϕi+1,j + ϕi−1,j + ϕi,j+1 + ϕi,j−1) .

Po želji v desno stran sproti vstavljamo že izračunane nove komponente. Lahko pa celo dobljeni novi približek ϕi1 poimenujemo kot "provizoričnega" in iz njega izračunamo "pravega",

ϕi2 = ϕi0 + ω(ϕi1ϕi0), 1 ≤ ω < 2 .

Kompleksna difuzijska enačba

Kompleksna difuzijska enačba (za spreminjanje amplitude stanja)

i

Ψ

t

=

2Ψ

x2

+ V(x)Ψ = HΨ

Implicitna metoda.

ima formalno rešitev Ψ(x,t) = eiHtΨ(x,0). Eksponentni operator aproksimiramo za kratek časovni korak z operatorjem "naprej-nazaj" (1 − 1/2 iHΔt) / (1 + 1/2 iHΔt), ki ohranja normiranost amplitude, in pridelamo enačbo v obliki

aΨ1i−1 + biΨ1i + aΨ1i+1 = F0i .

Količine a, bi in Fi so poznane funkcije začetne amplitude in potenciala. Taka enačba velja za vsako točko iz diskretnega nabora. Matrično zapisani sistem enačb je tridiagonalen in ga rešujemo z ustrezno eliminacijo.

Kompleksna potencialna enačba

Kompleksna potencialna enačba (za lastne energije in lastne amplitude stanja)

2ψ

dx2

= [EV(x)] ψ

Strelska metoda.

se zapiše v diskretnih točkah kot

ψi+1 = −ψi−1 + 2ψi  + F(E,Vi)ψi.

Izberemo primerno vrednost energije E. Ker mora veljavna amplituda izginiti pri prodiranju v visok potencial, postavimo začetno vrednost ψ0 na nič in naslednjo vrednost ψ1 na poljubno majhno število. Nato izračunamo vrednost ψ2 iz obeh predhodnih. Tako korakamo do desnega roba. Dobljeno amplitudo normiramo. S tem se ustrezno prilagodijo vse strmine, tudi tista, ki smo jo izbrali na levem robu. Če je sedaj desna robna amplituda enaka nič, je bila izbrana energija kar prava. Če pa ne, poskusimo znova z drugo energijo. Ko smo našli dve energiji, pri katerih amplituda spremeni predznak na desnem robu, poiščemo boljše približke z razpolavljanjem tega intervala.

Spektralna analiza

Vsako periodično funkcijo lahko zapišemo kot vsoto harmoničnih funkcij. Naj bo periodična funkcija f(t) poznana v ekvidistantnih točkah kΔt, k = 0, 1, 2, 3 … N − 1 preko celotne periode. Poznamo torej vrednosti fk. Funkcijo tedaj zapišemo kot superpozicijo:

fk =

1

N

N − 1

n=0

n e2πink/N .

Harmonične koeficiente izračunamo takole:

n =

N−1

k=0

  fk e−2πikn/N .

Če časovna funkcija ni periodična, moramo v interval na obeh koncih zajeti točke, ki so zanjo "tipične", karkoli pač to že pomeni. Ugodno je, če zna računalnik delati s kompleksnimi števili, sicer moramo za to poskrbeti sami z uporabo realnih dvojic.

Slučajna števila

Modeliranje slučajnih pojavov v naravi zahteva generacijo slučajnih števil s takšno ali drugačno porazdelitvijo.

Enakomerna porazdelitev.

Naravna števila z bolj ali manj enakomerno slučajno porazdelitvijo po intervalu [0, m−1] dobimo z iteracijo

xn+1 = (axn + b) mod m

iz primernega začetnega "semena" x0. Števila xn so tem bolj "slučajna", čim "boljša" naravna števila a, b in m uporabljamo za njihovo tvorjenje. Dobra je naslednja izbira:

a = 75 = 16 807
b = 0
m = 231 − 1 = 2 147 483 647 .

Če hočemo pridobivati realna slučajna števila na intervalu [0, 1], vsako dobljeno število delimo z modulom: rn = xn/m. Računalnik mora znati računati z dovolj velikimi naravnimi števili, namreč do am + b, sicer je potrebno algoritem prilagoditi.

Normalna porazdelitev.

Dvojica slučajnih števil x1 in x2 iz enakomerne porazdelitve po enotnem intervalu [0, 1] določa dvojico števil z1 in z2 iz normalne porazdelitve N(0, 1) takole:

y1 = 2x1 − 1
y2 = 2x2 − 1
r = y12 + y22
w = √(−2 ln r / r)
z1 = y1w
z2 = y2w .

Izračunani r mora biti manjši od 1, sicer postopek prekinemo in začnemo z novo dvojico začetnih števil. Normalno porazdelitev Z ∈ N(μ,σ) izračunamo iz standardne z ∈ N(0,1) s transformacijo Z = μ + σz.

M. Divjak