20

Numerika

Računalniki – Koreni enačb – Sistem linearnih enačb – Odvajanje – Integriranje – Spektralna analiza  – Enačba rasti – Enačba gibanja – Advekcijska enačba – Valovna enačba – Difuzijska enačba – Potencialna enačba – Amplitudna enačba

20.1 Računalniki

Matematika kot računsko orodje znanosti se ukvarja s števili, funkcijami in enačbami. V principu lahko vse to delamo s svinčnikom na papirju. Z izumom računalnika pa se vse spremeni. Današnji računalniki opravijo v sekundi toliko osnovnih računskih operacij, kolikor bi jih človek na papirju v milijon letih. Računi, ki so bili do sedaj preobsežni, postanejo praktično izvršjivi. Poglejmo, kako lahko računalnik uporabimo za reševanje tipičnih matematičnih problemov!

Slika 20.1 Osebni računalnik. Z njim komuniciramo preko tipkovnice in katodnega zaslona. (Anon)

20.2 Koreni 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.

Bisekcija

Z grobim tabeliranjem najdemo dve vrednosti a in b, pri katerih ima funkcija nasprotna predznaka: ničla (koren) leži tedaj nekje na intervalu [a,b]. Izračunamo osrednjo točko

in funkcijsko vrednost f(c) v njej. Odvisno od predznaka funkcije v tej točki je novi interval [a,c] ali [c,b]. Nadaljujemo, dokler ne skrčimo intervala na željeno majhnost. To je reševanje enačbe z bisekcijo in je zmeraj uspešno.

Regula falsi

Namesto da iščemo središčno točko intervala, je bolje iskati točko, kjer premica skozi obe krajiščni točki seka abscisno os. Tej premici rečemo sekanta. Podobna trikotnika z vrhoma pri presečišču povesta f(b)/(b − c) = f(a)/(a − c), torej

Nato izberemo pravega izmed obeh podintervalov ter ponovimo postopek. To je reševanje enačbe z metodo regula falsi. Potrebnih je manj korakov kot pri razpolavljanju.

Navadna iteracija

Morda lahko enačbo f(x) = 0 izrazimo kot x = g(x)? Potem vstavimo v desno stran primeren približek x₀ in izračunamo levo stran – novi približek x₁:

Tako nadaljujemo in upamo, da se bodo zaporedni približki vse bolj stiskali – konvergirali – k neki mejni vrednosti. Pravimo, da enačbo rešujemo z navadno iteracijo.

Kdaj pride do konvergence? Naj bo α iskani koren. Razvoj v vrsto pove g(x) = g(α) + (x − α)g'(α) + … Ker g(α) = α, velja g(x) − α ≈ g'(α)(x − α). Ker x_n+1 = g(x_n), sledi (x_n+1 − α) ≈ g'(α)(x_n − α). Drugače rečeno: razlika med približkom in korenom se pri vsaki iteraciji pomnoži približno z g'(α). Zato pride do konvergence, če |g'(α| < 1. Funkcija g(x) torej ne sme v okolici korena naraščati ali padati prestrmo.

20.3 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⁻¹ · b. Numerično jo določimo takole (Gauss).

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. To je metoda eliminacije.

Navadna 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 D = diag A (z ničlami drugod), druga pa preostanek R = A − D (z ničlami po diagonali). Tako dobimo enačbo D · x = b − R · x. V desno stran vstavimo primeren približek x⁰ in izračunamo levi vektor – z diagonalnimi koeficienti pomnoženi novi približek x¹. Eksplicitno velja:

Po želji lahko v desno stran sproti vstavljamo že izračunane nove komponente. Metoda konvergira, če je vsak diagonalni element matrike A po absolutni vrednosti večji od vsote absolutnih vrednosti preostalih elementov v vrstici. (Tako pravijo tisti, ki se na to spoznajo. Mi se dokazu hvaležno odrekamo.)

Prekomerna relaksacija

Ko iz približka x_i⁰ izračunamo novi približek x_i¹, je ta od predhodnika bolj ali manj različen. Naravno se zdi predpostaviti, da je "pravi" naslednji približek x_i² odvisen od prejšnjega približka x_i⁰ in razlike približkov x_i¹ − x_i⁰. Dobljeni približek x_i¹ torej vzamemo za "provizoričnega" in iz njega izračunamo "pravega":

To je prekomerna relaksacija. Če izberemo ω = 1, se enačba poenostavi v navadno iteracijo, kakor tudi mora biti. Parameter ω izberemo s poskušanjem tako, da je konvergenca čim hitrejša. Tipično je nekaj večji od 1. Pri vrednostih preko 2 pa rado pride do divergiranja.

20.4 Odvajanje

Kadar kakšne funkcije ne moremo ali nočemo odvajati algebraično, storimo to numerično, in sicer v vsaki točki, ki nas zanima.

Prvi odvod, enostranska shema

Okrog točke x, kjer iščemo odvod, razvijemo funkcijo u(x) v potenčno vrsto u(x + h) = u(x) + hu'(x) + …, iz nje izračunamo prvi odvod u' in zanemarimo vse višje člene:

To je napredna shema za izračun odvoda. Če bi razvili u(x − h) = u(x) − hu'(x) … , bi pa pridelali odzadnjo shemo

Napaka metode, ki jo pri razvoju obakrat zagrešimo, je sorazmerna s h: |u'_true − u'| ∝ h. Čim manjši h izberemo, tem manjša je ta napaka. Žal pa pri dejanskem računanju na končno število decimalnih mest h ne sme biti premajhen. Ker je namreč vsak člen v števcu obremenjen z zaokrožitveno napako, je relativna napaka njune razlike tem večja, čim manjši je h.

Prvi odvod, centralna shema

Ugodno bi bilo, ko bi bila napaka metode sorazmerna z višjo, ne zgolj s prvo potenco h. Očitno bo treba v razvoju funkcije v potenčno vrsto upoštevati več členov. Tako razvijemo u(x + h) in u(x − h), drugo enačbo odštejemo od prve in iz dobljene razlike izračunamo prvi odvod u', pri čemer zanemarimo vse višje člene:

To je centralna shema za izračun odvoda. Napaka metode je zdaj sorazmerna s h².

Drugi odvod, centralna shema

Po zgledu približkov za prvi odvod izračunajmo še približek za drugi odvod. Razvijemo u(x + h) in u(x − h), obe enačbi seštejemo, iz dobljene vsote izračunamo drugi odvod u" in zanemarimo višje člene, pa dobimo:

To je centralna shema za izračun drugega odvoda. Napaka metode je sorazmerna s h².

20.5 Integriranje

Ko odpovejo vsi algebraični prijemi za izračun integrala, ni druge, kot da se zatečemo k numeričnim metodam. Abscisno os razdelimo na primerne velike zaporedne intervale, nekako izračunamo delne ploščine nad vsakim in jih nato seštejemo.

Pravokotna shema

Prva misel je, da kos krivulje nad obdelovanim abscisnim intervalom [x, x + h] aproksimiramo s konstanto; tako pridelamo ploskev v obliki pravokotnika. Funkcijo u(x) torej razvijemo v potenčno vrsto okrog točke x in to do prvega člena; tako dobimo pravokotno shemo

Napaka izračunanega integrala je sorazmerna s h². Če hočemo integrirati na velikem intervalu [a, b], ga razdelimo na majhne intervale (b − a)/N = h. Levi robovi teh intervalov imajo koordinate x_i = a + ih . Celotni integral je potem vsota delnih integralov

Lokalne napake se pri seštevanju preko N = [b − a]/h podintervalov akumulirajo in celotna napaka N · h² postane sorazmerna z |b − a|h.

Trapezna shema

Druga misel je, da kos krivulje nad obdelovanim abscisnim intervalom aproksimiramo s sekanto; tako pridelamo ploskev v obliki trapeza. Funkcijo u(x) torej razvijemo v potenčno vrsto okrog točke x do linearnega člena, izrazimo prvi odvod z enostransko desno shemo, zanemarimo višje člene in dobimo – z napako, sorazmerno s h³ – trapezno shemo

Integral po večjem območju [a, b] je vsota integralov po podobmočjih (b − a)/N = h. Označimo x_i = a + ih, pa dobimo

Lokalne napake se pri seštevanju akumulirajo in celotna napaka postane sorazmerna z |b − a|h².

Parabolična shema

Še bolj natančna je aproksimacija s parabolo. Funkcijo torej razvijemo v potenčno vrsto okrog točke x + h do kvadratnega člena, izrazimo prvi odvod s centralno shemo, drugi odvod s centralno shemo, zanemarimo višje člene in dobimo – z napako, ki je zdaj sorazmerna s h⁵ – parabolično shemo (Simpson)

Integral po večjem območju [a, b] je vsota integralov po podobmočjih (b − a)/2N = h, torej

Kumulativna napaka je sorazmerna z |b − a|h⁴. Od vseh naštetih shem je torej parabolična daleč najbolj natančna in je zato priporočljiva za uporabo.

20.6 Spektralna analiza

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

Harmonične koeficiente izračunamo takole:

Č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.

Vzorčenje časovne funkcije z odseki Δt ne more zajeti harmonikov, ki imajo periodo krajšo od 2Δt, ampak jih prepozna kot harmonike z daljšimi periodami. Zato tudi ni verodostojno računati njihovih amplitud. To je razlog, zakaj pri diskretni transformaciji računamo največ toliko harmonikov, kot imamo na voljo časovnih točk.

Filtriranje šuma

Vremenoslovci merijo temperaturo zraka vsak dan ob treh časih: zjutraj, opoldne in zvečer. Tako dobijo časovni niz preko mnogih let. Ta niz kaže počasno nihanje preko enega leta (sezonske spremembe), na katerega je naloženo hitro nihanje preko enega dne (dnevne spremembe). Kako bi iz časovnega niza "odstranili" dnevne spremembe, da ne bi "kvarili" sezonskih sprememb oziroma da bi bile te bolj vidne? — Časovnemu nizu določimo spekter. — Dobljeni spekter pomnožimo s tako funkcijo, da nam frekvence izven izbranega frekvenčnega pasu zavzamejo vrednost 0, znotraj pa ostanejo nespremenjene. — Iz popravljenega spektra izračunamo popravljeni časovni niz, ki ne vsebuje več motečih frekvenc. Rečemo, da smo niz filtrirali. Opisana metoda je zelo priročna za filtriranje vsakršnih nizov, v katerih je signal pomešan s šumom.

20.7 Enačba rasti

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

Tangentna metoda

rešujemo najbolj preprosto takole. Odvod dm/dt izrazimo z napredno shemo (20.6) in dobimo

To je tangentna metoda (Euler). Če poznamo vrednost m(t) ob času t, lahko iz (20.19) izračunamo, kakšna je vrednost m(t + dt) ob malo kasnejšem času t + dt. Nato postopek ponavljamo po korakih dt. Za zagon potrebujemo začetni pogoj m(t₀) = m₀.

Metoda ima – zaradi uporabljene sheme odvoda – pri enem koraku lokalno napako |m_true − m| ∝ (dt)². Čim manjši korak dt uporabimo, tem bolj natančno je določena m(t + dt). Premajhnega koraka pa spet nima smisla vzeti, ker napako (dt)² prevpije nenatančnost računanja na N decimalnih mest, 10^−N. Za najmanjši še smiselni časovni korak zato velja dt > 10^−N/2.

Z naraščanjem števila korakov se lokalne napake akumulirajo: kumulativna napaka po n korakih je zato sorazmerna z n · (dt)² = |t − t₀|/dt · (dt)² = |t − t₀|dt. Tetivna metoda je zato malo natančna. Uporabna je predvsem za računanje ne predaleč od začetne točke.

Dvotangentna metoda

Poskusimo najti boljšo shemo! Funkcijo m(t + dt) razvijemo v vrsto do kvadratnega člena. V tem razvoju upoštevamo dm/dt = f(m,t) in d²m/dt² = df(m,t)/dt. Nato aproksimiramo df(m,t)/dt = [f(m(t + dt),t + dt) −f(m,t)]/dt in znotraj tega m(t + dt) = m(t) + dt f(m,t). Ko vse skupaj zložimo, dobimo

Namesto tangente dm/dt v točki t torej uporablja metoda povprečno vrednost tangent v točkah t in t+dt. Da bi določili vrednost tangente v točki t+dt, bi pravzaprav že morali poznati vrednost m(t+dt) v tej točki. Ker tega ne vemo, aproksimiramo vrednost m(t+dt) s formulo (20.6). Za praktično računanje so primerne naslednje oznake:

To je dvotangentna metoda. Njena napaka pri enem koraku je sorazmerna z (dt)³ in preko daljšega intervala sorazmerna s |t − t₀|(dt)².

20.8 Enačba gibanja

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

je ekvivalentna sistemu dveh sklopljenih enačb

Dvotangentna metoda

To sta dve enačbi prvega reda in rešujemo ju skupaj, korak za korakom, po dvotangentni metodi (20.21): najprej izračunamo v(t + dt), nato pa še s(t + dt). Lego in hitrost torej računamo v istih časovnih točkah.

Preskočna metoda

Računanje lege in hitrosti v istih časovnih točkah gotovo ni najbolje. Če se omejimo na primer, ko sila f ni odvisna od hitrosti v, sta lega in hitrost lepo "prepleteni" v času: s(t + dt) = s(t) + v(t + dt/2) dt in v(t + dt/2) = v(t − dt/2) + f(t,s(t))dt. Obe količini lahko zato računamo v medsebojno zamaknjenih časovnih točkah. Ob času t naj bo začetna lega s₀ in ob času t + dt/2 "začetna" hitrost v_1/2. Nove vrednosti dt kasneje so:

Če "začetne" hitrosti ne poznamo, jo določimo iz prave začetne hitrosti v₀ s posebnim korakom po tangentni metodi v_1/2 = v₀ + f₀dt/2.

20.9 Advekcijska enačba

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

preučujmo na intervalu [0,l]. Interval opremimo s točkami v razmikih dx. V teh točkah si mislimo vrednosti Q_iⁿ ob času ndt. Iz njih moramo izračunati vrednosti Q_iⁿ⁺¹ ob naslednjem času (n+1)dt.

Protitočna shema

V vsaki točki aproksimiramo časovni odvod z diferenco naprej v času in prostorski odvod z diferenco proti toku: ∂Q/∂t = (Qⁿ⁺¹ − Qⁿ)/dt in ∂Q/∂x = (Q_i − Q_i−1)/dx. Nove vrednosti so potem podane eksplicitno takole (za c > 0):

Računati začnemo ob času n = 0, ko so točke opremljene z začetnim stanjem. Za vsako notranjo točko izračunamo prihodnjo vrednost. Robni točki izračunamo posebej, v skladu z robnimi pogoji: postavimo ju na predpisano vrednost ali na vrednost prve notranje točke. Potem nadaljujemo z naslednjim korakom v času, dokler je pač treba.

Kriterij stabilnosti

Kolikšna intervala dx in dt naj izberemo? Izkušnje kažejo, da pri uporabljenem r kakšna točkovna vrednost sčasoma podivja v neskončnost. Ko primerno zmanjšamo r, pa se to ne zgodi. Da bo shema stabilna, je očitno potreben naslednji pogoj: max_i |Q_iⁿ⁺¹| ≤ max_i |Q_iⁿ|. Analitična rešitev Q(x,t) advekcijske enačbe je vsak izraz oblike exp (iωt) exp(ikx), pa tudi linearna kombinacija takih členov. (Paziti moramo na razliko med kompleksno enoto i in točko i.) Naj bo torej Q_iⁿ = exp(iωndt) exp(ikidx) takšna elementarna rešitev v točkah (i,n). Potem max_i|Q_iⁿ| = |exp(iωndt)| in kriterij stabilnosti se zapiše kot

Stabilnost protitočne sheme (20.26) torej raziščemo tako, da vanjo vstavimo Q_iⁿ⁺¹ = exp(iω(n+1)dt in Q_iⁿ = exp(iωndt) ter izračunamo njun količnik. Ta vsebuje parameter r in s kriterijem stabilnosti je slednji tudi določen. Tako izračunamo λ = 1 − r(1 − cos kdx) − ir sin kdx) in iz kriterija |λ| ≤ 1 sledi

Torej moramo uporabiti tako kratek dt, da se lokalne motnje premaknejo za manj kot dx. Če hočemo zmanjšati prostorski interval za dvakrat, moramo zmanjšati časovni interval za dvakrat, torej povečati računsko delo za faktor štiri.

20.10 Valovna enačba

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

se zapiše kot sistem dveh advekcijskih enačb

Nazorno sta to enačbi za gladinske valove v plitvi vodi: h je višina vala nad/pod neperturbirano gladino (normirana na njeno globino), v pa horizontalna hitrost vode. Gradient hitrosti torej povzroča spremembo višine vala, gradient višine vala pa povzroča spremembo hitrosti.

Preskočna shema

Zaradi prepletenosti obeh spremenljivk h in v in njunih gradientov je smiselno računanje v dveh naborih točk. Naj bosta torej začetni polji v dveh naborih točk v⁰_i in h⁰_i+1/2; nabora točk sta medsebojno zamaknjena za interval dx / 2. V časovnem intervalu dt se hitrostno polje spremeni v

nakar se spremeni še višina polja v

To je preskočna shema. Na robovih je treba upoštevati primerne robne pogoje. Zaradi medsebojne zamaknjenosti točk so odvodi efektivno centralni in metoda je drugega reda natančnosti, to je, njena kumulativna napaka je sorazmerna s |t−t₀|(dt)².

Kriterij stabilnosti

Da bo rešitev dveh advekcijskih enačb stabilna, moramo uporabljati, kot smo že spoznali, dovolj kratke časovne korake:

20.11 Difuzijska enačba

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

Napredno-centralna shema

rešujemo v točkah i z diferencami naprej v času in centralno v prostoru:

Računati začnemo ob času n = 0, ko so točke opremljene z začetnim stanjem. Za vsako notranjo točko izračunamo prihodnjo vrednost. Robni točki izračunamo posebej, v skladu z robnimi pogoji: postavimo ju na predpisano vrednost ali na vrednost prve notranje točke. Potem nadaljujemo z naslednjim korakom v času, dokler je pač treba.

Kriterij stabilnosti

Stabilnost sheme določimo tako kot pri advekcijski (in valovni) enačbi. V shemo (20.35) vstavimo Q_iⁿ⁺¹ = exp(iω(n+1)dt in Q_iⁿ = exp(iωndt) ter izračunamo njun količnik. Ta vsebuje parameter r in s kriterijem stabilnosti je slednji tudi določen. Tako izračunamo λ = 1 − 4r sin² (k/2) in iz kriterija |λ| ≤ 1 sledi

Če hočemo zmanjšati prostorski interval za dvakrat, moramo zmanjšati časovni interval za štirikrat, torej povečati računsko delo za faktor osem. V treh dimenzijah ravnamo podobno. Prostor razdelimo na kocke z robom dl in računamo Q_ijk v njihovih ogliščih. Stabilnost sedaj zahteva r = Ddt/dl² ≤ 1/6.

20.12 Potencialna enačba

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

Prekomerna relaksacija

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

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

Robni pogoji

Pri vsakem časovnem koraku morajo biti podani vsi robni pogoji s predpisanimi vrednostmi. Če je kakšen robni pogoj podan z odvodom, na primer s ϕ_x(i=1) = A, izračunamo robno vrednost ϕ(1) iz aproksimacije A = [ϕ(2) − ϕ(1)]/dx. Posebej za ϕ_x(i=1) = 0 velja ϕ(1) = ϕ(2).

Področje, na katerem želimo rešiti potencialno enačbo, tudi ni nujno pravokotnik. Če je krog ali krogla, si pomagamo tako, da zapišemo ∇² ϕ = 0 v cilindričnih ali krogelnih koordinatah in primerno diskretiziramo odvode. Če pa je področje "nepravilne" oblike, potem vozlišča mreže ne padejo točno na robove območja. Za vsako robno točko mreže moramo potem določiti vrednost z interpolacijo iz notranjih točk. V podrobnosti se ne bomo spuščali.

20.13 Amplitudna enačba

Gibanje kvantnega delca v potencialnem polju opisuje razširjena amplitudna enačba

Energije E ne poznamo. V diskretnih točkah zapišemo enačbo kot

Strelska metoda.

Nato izberemo primerno vrednost energije E. Ker mora veljavna amplituda izginiti pri prodiranju v visok potencial, postavimo začetno vrednost ψ₀ na nič in naslednjo vrednost ψ₁ na poljubno majhno število. Nato izračunamo vrednost ψ₂ 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. Podobno poiščemo tudi druge energije in lastne funkcije. Metoda je uporabna tudi za osnosimetrične in centralne potenciale V(r), le ∇²ψ moramo zapisati v ustreznih koordinatah. □