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.
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 + θ(u2 − u1) . |
Parabolična interpolacija
Bolj natančna je aproksimacija s parabolo skozi tri zaporedne ekvidistantne točke:
| u = u1 + θ(u2 − u1) + (1/2)θ(θ − 1)(u3 − 2u2 + u1) . |
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 = a − f(a) | b − a f(b) − f(a) |
. |
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:
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 = b − R · 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 − |
∑ j ≠ i |
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 + ω(xi1 − xi0), 1 ≤ ω < 2 . |
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φ ≡ | Aqq − App 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 r ≠ p, r ≠ q):
| A'pq = 0 |
| A'pp = App − tApq |
| A'qq = Aqq − tApq |
| A'rp = Arp − s(Arq + τArp) |
| A'rq = Arp − s(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.
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(x − h) 2h |
d2u dx2 |
≈ | u(x + h) − 2u(x) + u(x − h) 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č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 (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.
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 = Qi − r(Qi+1 − Qi) |
| 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 + Qi − Qi+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 = bi − mci−1 |
| di = di − mdi−1 , |
nato pa za člene od n−1 do 1 po vrsti nazaj
| xi = (di − cixi+1) / bi . |
Shema je stabilna za poljuben časovni korak.
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)Qi1 − rQi+11 = rQi−1 + (1 − 2r)Qi + rQi+1 |
| r = Ddt / 2dx2 . |
Shema je stabilna za vsak časovni korak.
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 = E0i − c2 (dt / dx) (B0i+1/2 − B0i−1/2) , |
nakar se spremeni še magnetno polje v
| B1i+1/2 = B0i+1/2 − (dt / dx) (E1i+1 − E1i) . |
Na robovih je treba upoštevati primerne robne pogoje. Da bo rešitev stabilna, moramo uporabljati dovolj kratke časovne korake: cdt / dx ≤ 1.
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 (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 (za lastne energije in lastne amplitude stanja)
∂2ψ dx2 |
= [E − V(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.
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 |
B̂n e2πink/N . |
Harmonične koeficiente izračunamo takole:
| B̂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.
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.