Numerično odvajanje

Uvod¶

Vsako elementarno funkcijo lahko analitično odvajamo. Definicija odvoda je:

Neposredna uporaba zgornje enačbe vodi v odštevanje zelo podobnih funkcijskih vrednostih ($f(x+\Delta x)$, $f(x)$), obremenjenih z zaokrožitveno napako, ki jih delimo z majhno vrednostjo $\Delta x$; posledično ima odvod bistveno manj signifikantnih števk kakor pa funkcijske vrednosti. Numeričnemu odvajanju se izognemo, če imamo to možnost; je pa v nekaterih primerih (npr. reševanje diferencialnih enačb) nepogrešljivo orodje!

Pri numeričnem odvajanju imamo dva, v principu različna, pristopa:

1. najprej izvedemo interpolacijo/aproksimacijo, nato pa na podlagi znanih interpolacijskih/aproksimacijskih funkcij izračunamo odvod (o tej temi smo že govorili pri interpolaciji oz. aproksimaciji) in

2. računanje odvoda neposredno iz vrednosti iz tabele.

V okviru tega poglavja se bomo seznanili s tem, kako numerično izračunamo odvod funkcije $f(x)$; pri tem so vrednosti funkcije $f(x)$ podane tabelarično (pari $x_i$, $y_i$), kakor je prikazano na sliki: Najprej se bomo osredotočili na ekvidistantno, s korakom $h$, razporejene vrednosti $x_i$; vrednosti funkcije pa bodo $y_i=f(x_i)$.

Glede na zgornjo definicijo odvoda, bi prvi odvod (za mesto $i$) lahko zapisali:

kjer je $h=x_{i+1}-x_{i}$. S preoblikovanjem enačbe:

lahko tudi rečemo, da za prvi odvod funkcije na mestu $i$, utežimo funkcijsko vrednost pri $i$ z $-1/h$ in funkcijsko vrednost pri $i+1$ z $+1/h$.

V nadaljevanju si bomo pogledali teoretično ozadje kako določimo ustrezne uteži za različne stopnje odvodov, katere možnosti pri tem imamo in kako to vpliva na red natančnosti.

Aproksimacija prvega odvoda po metodi končnih razlik¶

Za uvod si oglejmo spodnji video:

Odvod $f'(x)$ lahko aproksimiramo na podlagi razvoja Taylorjeve vrste. To metodo imenujemo metoda končnih razlik ali tudi diferenčna metoda.

Razvijmo Taylorjevo vrsto naprej (naprej, zaradi člena $+h$):

Člen $\mathcal{O}\left(h^{2}\right)$ označuje napako drugega reda. Če iz enačbe izrazimo prvi odvod:

Ugotovimo, da lahko ocenimo prvi odvod v točki $x_i$ (to je: $f_o'(x_i)$) na podlagi dveh zaporednih funkcijskih vrednosti:

in pri tem naredimo napako metode, ki je prvega reda $\mathcal{O}\left(h^{1}\right)$.

Uporabili smo $y_i=f(x_i)$ (glejte sliko zgoraj).

Napaka je:

kjer je $\xi$ neznana vrednost na intervalu $[x_i, x_{i+1}]$ in smo zanemarili višje člene.

Velja torej izraz:

Sedaj si poglejmo, kako pridemo do istega rezultata s strojno izpeljavo; najprej uvozimo sympy:

Definirajmo simbole:

Nato nadaljujemo z razvojem Taylorjeve vrste naprej (angl. forward Taylor series):

Člen $\mathcal{O}\left(h^{2}\right)$ vsebuje člene drugega in višjega reda. V zgornji enačbi je uporabljena začasna spremenljivko za odvajanje $\xi_1$; izvedmo odvajanje in vstavimo $\xi_1=x$:

Zapišemo enačbo:

Rešimo jo za prvi odvod $f'(x)$:

V kolikor odvoda drugega in višjih redov ne upoštevamo, smo naredili torej napako:

Napaka $\mathcal{O}\left(h^{1}\right)$ je torej prvega reda in če ta člen zanemarimo, naredimo napako metode in dobimo oceno odvoda:

Ugotovimo, da gre za isti izraz, kakor smo ga izpeljali zgoraj, torej je:

Uteži torej so:

Centralna diferenčna shema¶

Odvod $f'(x)$¶

Najprej si poglejmo razvoj Taylorjeve vrste nazaj (angl. backward Taylor series):

Ugotovimo, da se pri razliki vrste naprej in nazaj odštevajo členi sodega reda; definirajmo:

Izvedemo sledeče korake:

1. Taylorjevo vrsto nazaj odštejemo od vrste naprej, sodi odvodi se odštejejo,

2. rešimo enačbo za prvi odvod,

3. določimo napako metode,

4. določimo oceno odvoda.

Izvedimo zgornje korake:

Ocena 1. odvoda torej je:

Ali:

Uteži torej so:

Napaka metode pa je torej drugega reda:

Zgled: $\exp(-x)$¶

Poglejmo si zgled eksponentne funkcije $f(x)=\exp(-x)$ in za točko $x=1,0$ izračunajmo prvi odvod $f'(x)=-\exp(-x)$ pri koraku $h_0=1$ in $h_1=0,1$.

Najprej pripravimo tabelo numeričnih vrednosti in točen rezultat:

Potem uporabimo shemo naprej:

Izračunajmo napako pri $x=1,0$:

Potrdimo lahko, da je napaka pri koraku $h/10$ res približno 1/10 tiste pri koraku $h$.

Poglejmo sedaj še napako za centralno diferenčno shemo, ki je drugega reda:

Analizirajmo napako:

Potrdimo lahko, da je napaka pri koraku $h/10$ res približno 1/100 tiste pri koraku $h$.

Odvod $f''(x)$¶

Če Taylorjevo vrsto naprej in nazaj seštejemo, se odštejejo lihi odvodi:

Določimo drugi odvod:

Ocena drugega odvoda je:

Ali:

Uteži torej so:

Napaka metode pa je ponovno drugega reda:

Odvod $f'''(x)$¶

Če želimo določiti tretji odvod, moramo Taylorjevo vrsto razviti do stopnje 5:

Uporaba 1. odvoda, ki smo ga izpeljali zgoraj, nam ne bi koristila, saj je red napake $\mathcal{O}\left(h^{2}\right)$, kar pomeni, da bi v zgornji enačbi pri deljenju s $h^3$ dobili $\mathcal{O}\left(h^{-1}\right)$.

Uporabimo trik: ponovimo razvoj, vendar na podlagi dodatnih točk, ki sta od $x$ oddaljeni za $2h$ in $-2h$:

Sedaj imamo dve enačbi in dve neznanki; sistem bomo rešili po korakih:

1. enačbo eq_h rešimo za prvi odvod,

2. enačbo eq_2h rešimo za prvi odvod,

3. enačimo rezultata prvih dveh korakov in rešimo za tretji odvod,

4. določimo napako metode,

5. določimo oceno odvoda.

Izvedimo navedene korake:

Ocena 3. odvoda je:

Ali:

Uteži torej so:

Potrdimo, da je napaka metode drugega reda:

Odvod $f^{(4)}(x)$¶

Ponovimo podoben postopek kot za 3. odvod, vendar za 4. odvod seštevamo Taylorjevo vrsto (do stopnje 6) naprej in nazaj:

Pripravimo dodatno enačbo na podlagi točk, ki sta od $x$ oddaljeni za $2h$ in $-2h$:

Iz dveh enačb določimo 4. odvod:

Ocena 4. odvoda je:

Ali:

Uteži torej so:

Potrdimo, da je napaka metode drugega reda:

Povzetek centralne diferenčne sheme¶

Zgoraj smo izpeljali prve štiri odvode z napako metode 2. reda. Bistvo zgornjih izpeljav je, da nam dajo uteži, s katerimi moramo množiti funkcijske vrednosti, da izračunamo približek določenega odvoda. Iz tega razloga bomo tukaj te uteži zbrali. Če ste neučakani, lahko skočite na tabelo spodaj. Z branjem nadaljujte, če pa želite spoznati, kako predhodno izpeljane izraze strojno uredimo.

Najprej zberimo vse ocene odvodov v seznam:

Na razpolago imamo 5 funkcijskih vrednosti (pri legah $x-2h, x-h, x, x+h, x+2h$), ki jih damo v seznam:

Utež prvega odvoda za funkcijsko vrednosti $f(x-h)$ izračunamo:

Sedaj posplošimo in izračunajmo uteži za vse funkcijske vrednosti in za vse ocene odvodov:

Zgornje povzetke lahko tudi zapišemo v tabelarični obliki:

Prikazana centralna diferenčna shema ima napako 2. reda $\mathcal{O}(h^{2})$.

Opomba: v kolikor vas zanima posplošitev diferenčne metode, prosim glejte paket findiff.

Izboljšan približek - Richardsonova ekstrapolacija¶

Če je točen odvod je izračunan kot:

kjer je $f_o'(x_i)$ numerično izračunan odvod v točki $x_i$ in $e$ ocena napake.

Za metodo reda točnosti $n$: $\mathcal{O}(h^{n})$ pri koraku $h$ velja:

kjer je $K$ neznana konstanta.

Če korak razpolovimo in predpostavimo, da se $K$ ne spremeni, velja:

Iz obeh enačb izločimo konstanto $K$ in določimo izboljšan približek:

Zgled¶

Poglejmo si zgled $f(x)=\sin(x)$ (analitični odvod je: $f'(x)=\cos(x)$):

Pri koraku $h$ imamo funkcijske vrednosti definirane pri:

Numerični odvod pri $x=\pi$:

in koraku $h$ je

Izračun pri koraku $2h$:

Izračunajmo izboljšano oceno za $x=\pi$:

Vidimo, da je izboljšana ocena najbližje teoretični vrednosti $\cos(\pi)=-1$.

Necentralna diferenčna shema¶

Centralna diferenčna shema, ki smo jo spoznali zgoraj, je zelo uporabna in relativno natančna. Ker pa je ne moremo vedno uporabiti (recimo na začetku ali koncu tabele), si moramo pomagati z necentralnimi diferenčnimi shemami za računanje odvodov.

Poznamo:

• diferenčno shemo naprej, ki odvod točke aproksimira z vrednostmi funkcije v naslednjih točkah in

• diferenčno shemo nazaj, ki odvod točke aproksimira z vrednostmi v predhodnih točkah.

Izpeljave so podobne, kakor smo prikazali za centralno diferenčno shemo, zato jih tukaj ne bomo obravnavali in bomo prikazali samo končni rezultat.

Diferenčna shema naprej¶

Diferenčna shema naprej z redom napake $\mathcal{O}(h^{1})$:

Diferenčna shema naprej z redom napake $\mathcal{O}(h^{2})$:

Diferenčna shema nazaj¶

Diferenčna shema nazaj z redom napake $\mathcal{O}(h^{1})$:

Diferenčna shema nazaj z redom napake $\mathcal{O}(h^{2})$:

Uporaba numpy.gradient¶

Za izračun numeričnih odvodov (centralna diferenčna shema 2. reda) lahko uporabimo tudi numpy.gradient() (dokumentacija):

gradient(f, *varargs, **kwargs)

kjer f predstavlja tabelo vrednosti (v obliki numeričnega polja) funkcije, katere odvod iščemo. f je lahko ene ali več dimenzij. Pozicijski parametri varargs definirajo razdaljo med vrednostmi argumenta funkcije f; privzeta vrednost je 1. Ta vrednost je lahko skalar, lahko pa tudi seznam vrednosti neodvisne spremenljivke (ali tudi kombinacija obojega). Gradientna metoda na robovih uporabi shemo naprej oziroma nazaj; parameter edge_order definira red sheme, ki se uporabi na robovih (izbiramo lahko med 1 ali 2, privzeta vrednost je 1).

Rezultat funkcije gradient je numerični seznam (ali seznam numeričnih seznamov) z izračunanimi odvodi.

Za podrobnosti glejte dokumentacijo.

Za odvajanje funkcije (in ne tabele vrednosti), glejte: scipy.differentiate.

Zgled¶

Pogledali si bomo zgled, kako uporabimo uteži, funkcijo gradient in posebnosti na robovih. Najprej pripravimo tabelo podatkov:

Uteži diferenčih shem:

Sedaj izvedemo odvod notranjih točk (prvi način je z izpeljevanjem seznamov, drugi je vektoriziran):

Na robovih uporabimo diferenčno shemo naprej oziroma nazaj:

Sestavimo rezultat:

Prikažemo rezultat skupaj z rezultatom funkcije np.gradient:

Uporaba numpy.diff¶

Dajmo tukaj predstaviti še pogosto uporabljeno funkcijo numpy.diff (dokumentacija) izračuna n-to diskretno razliko vzdolž podane osi. Prva razlika je podana kot $a[i+1] - a[i]$. To je osnovna operacija za numerično odvajanje (gre za odvod prve stopnje natančnosti), vendar moramo rezultat deliti s korakom $h$.

Pri uporabi np.diff se dolžina polja zmanjša za 1.

Zaokrožitvena napaka pri numeričnem odvajanju¶

Zgoraj smo se osredotočili na napako metode. Pri numeričnem odvajanju pa moramo biti zelo pozorni tudi na zaokrožitveno (ali tudi upodobitveno) napako! Poglejmo si prvi odvod (po centralni diferenčni shemi) zapisan z napako metode ($k\,h^2$) in zaokrožitveno napako $\varepsilon$:

V najslabšem primeru se zaokrožitvena napaka sešteje in je skupna napaka:

Ko je $h$ velik prevladuje napaka metode $k\,h^2$; ko pa je $h$ majhen, pa prevladuje zaokrožitvena napaka. Napaka ima minimum, ko velja:

Sledi:

Zgled¶

Spodaj si bomo pogledali primer, kjer bomo natančnost spreminjali v treh korakih:

1. float16 - 16-bitni zapis: predznak 1 bit, 5 bitov eksponent, 10 bitov mantisa

2. float32 - 32-bitni zapis: predznak 1 bit, 8 bitov eksponent, 23 bitov mantisa

3. float64 - 64-bitni zapis: predznak 1 bit, 11 bitov eksponent, 52 bitov mantisa (to je privzeta natančnost).

Za več o tipih v numpy glejte dokumentacijo

Določimo sedaj osnovno zaokrožitveno napako za posamezni tip:

Kot primer si poglejmo seštevanje: k številu 1. prištejemo polovico osnovne zaokrožitvene napake eps16 in pretvorimo v tip float16, ugotovimo, da je nova vrednost še vedno enaka vrednosti 1.:

Definirajmo najprej funkcijo $\exp(x)$, ki bo dala rezultat natančnosti, ki jo definira parameter dtype:

Definirajmo še funkcijo za analitično določljiv odvod (to bomo pozneje potrebovali za določitev relativne napake):

Podobno kakor zgoraj pri seštevanju, lahko tudi pri vrednosti funkcije ugotovimo, da sprememba vrednosti $x$, ki je manjša od $\epsilon$, vodi v isti rezultat:

Uporabimo sedaj centralno diferenčno shemo za prvi odvod. Pri tem naj bodo števila zapisana z natančnostjo dtype, s pomočjo točnega odvoda pa se izračuna še relativna napaka:

Poglejmo primer odvoda pri $x=1,0$ (Python funkcija vrne vrednost in relativno napako):

Definirajmo sedaj korak:

Izračunamo oceno odvodov za različne natančnosti zapisa (zaradi deljenja z 0 dobimo opozorilo):

Izrišemo različne tipe v odvisnosti od velikosti koraka $h$. Najprej uvozimo potrebne knjižnice:

Definirajmo sliko:

Prikažimo jo:

Pri relativno velikem koraku $h$ prevladuje napaka metode, pri majhnem koraku pa zaokrožitvena napaka; optimalni korak lahko ocenimo glede na:

V konkretnem primeru velja:

Sledi:

Izračunamo primeren korak za 16, 32 in 64-bitni zapis:

Najbolje se izkaže 64-bitni zapis, vendar pa tudi pri tem korak manjši od cca 1e-5 ni priporočen!

Aproksimacija odvoda s kompleksnim korakom¶

Aproksimacija odvoda s kompleksnim korakom je uporabna takrat, ko imamo definirano katerokoli analitično funkcijo, vendar pa nimamo na voljo analitičnega izraza za prvi odvod, za podrobnosti glejte (vir). Ideja izhaja iz razvoja funkcije $f(x)$ v Taylorjevo vrsto, vendar se izvede korak v smeri imaginarne osi $i\,h$ ($i=\sqrt{-1}$):

Če sedaj predpostavimo, da funkcija $f(x)$ realne vrednosti slika na realno os in da sta $x$ in $h$ realni vrednosti, potem lahko izpeljemo:

Iz prve enačbe (zgoraj) potem izpeljemo izraz za prvi odvod:

Rezultat je na nek način zelo presenetljiv. Zgoraj smo z eno koračno shemo naprej uspeli pridobiti natančnost $O\left(h^{1}\right)$, tukaj pa $O\left(h^{2}\right)$!

Pomembna prednost metode s kompleksnim korakom je, da nima težav z odštevanjem skoraj enakih števil (ang. cancellation error), ki pesti klasične diferenčne metode pri zelo majhnih korakih $h$. Zato lahko izberemo izjemno majhen $h$ (npr. 10^-200) in dobimo praktično točen rezultat (do strojne natančnosti).

Če pogledamo sedaj uporabo metoda na primeru od zgoraj $\exp(-x)$, spomnimo se najprej točnega rezultata pri vrednosti $x=1$:

Pri klicu s korakom $h=0.1$ dobimo točni dve (tri) števki (podobno kakor pri centralni diferenčni shemi), zmanjšanjem koraka na desetino pa bi število točnih števk približno podvojili.

Ker nimamo težav z odštevanjem skoraj enakih števil, lahko korak izjemno zmanjšamo in dobimo točen rezulta (do strojne natančnosti):

Poglejmo si še primer odvoda funkcije sin(x**2). Iz kode in slike spodaj vidimo bistveno hitrejšo konvergenco metode kompleksnega koraka. Zakaj metoda deluje tako dobro? Če vstavimo $x+i\,h$ v $x^2$, dobimo: $(x+i\,h)^2 = x^2 + 2i\,x\,h - h^2$

Ker je $h$ zelo majhen, je $h^2$ zanemarljiv v realnem delu, imaginarni del $2\,x\,h$ pa se linearno prenese skozi sinus. Na tak način dobimo približek teoretičnega odvoda $2\,x\,\cos(x^2)$. Metoda vidi odvod v imaginarnem delu brez odštevanja, zato ne pride do napake zaokroževanja.

Dodatno¶

Video predavanja na temo numeričnega odvajanja.

Vprašanja za vaje¶

Vprašanje 1: Podane so izmerjene vrednosti opravljene poti tekača s pri $N$ točkah v času. Numerično izračunajte potek hitrosti in pospeška (uporabite lahko poljubno metodo).

Opazujte in komentirajte kaj se dogaja, ko povečujete število točk $N$.

Vprašanje 2: Simbolno je podana funkcija $f(x)$ in njen odvod, $f'(x)$.

Potek funkcije in njenega odvoda opazujemo pri 100 točkah časa $t\in [0, 1]$.

Pripravi funkciji x_fun, dx_fun za numerično računanje. Izriši pripravljeni funkciji na intervalu $t\in [0, 1]$ z uporabo paketa matplotlib.

Vprašanje 3: Razdelite opazovan časovni interval $t \in [0,1]$ na 100 ekvidistantnih segmentov.

Pripravite vektor vrednosti $x(t)$ in odvod $x'(t)$ izračunajte še numerično, z uporabo funkcije numpy.gradient(x, h).

Vprašanje 4: Uporabite funkcijo np.gradient še za izračun drugega odvoda $x''(t)$. Opazujte vpliv parametra edge_order na natančnost odvodov na robovih intervala.

Vprašanje 5: Numeričnemu polju $x(t)$ pri 500 ekvidistantnih točkah na intervalu $t \in [0, 1]$ dodajte naključni šum $n(t)$ amplitude $A=0.1$ z uporabo funkcijo np.random.randn.

Z uporabo funkcije np.gradient numerično določite vrednost prvega in drugega odvoda šumnega signala.

Uteži diferenčnih shem s paketom findiff¶

Primer 6: Izračunajte vrednosti drugega odvoda $x''(t)$. Uporabite centralno diferenčno shemo drugega odvoda čez 5 točk (red napake 4).

Vprašanje 7: Na podlagi numeričnih podatkov $x(t)$ pri delitvi intervala $t \in [0,1]$ na 100 točk izračunajte vrednost prvega odvoda, $x'(t)$ po enačbi centralne diferenčne sheme preko treh točk. Ali lahko na ta način določite vrednost odvoda v točkah $t=0$ in $t=1$?

Vprašanje 8: Z uporabo diferenčne sheme naprej oz. nazaj izračunajte prvi odvod $x'(t)$ tudi v točkah na robu ($t=0$, $t=1$). Primerjajte dobljene vrednosti z rezultatom funkcije np.gradient, kjer nastavite argument edge_order=2.

(Namig: uporabite lahko funkcijo np.dot)

Izboljšanje približka odvoda z Richardsonovo ekstrapolacijo¶

Vprašanje 9: Izračunajte izboljšan približek odvoda funkcije:

pri vrednosti $t_i = 0.5$, tako, da uporabite zgornjo enačbo za dve različni delitvi, $h=0.05$ in $h/2 = 0.025$. Uporabite centralno diferenčno shemo za 1. odvod za red napake $n=2$.

Uporaba np.convolve¶

Vprašanje 10: Uporabite funkcijo np.convolve in izračunajte vrednosti drugega odvoda $x(t)$ pri delitvi intervala $t \in [0, 1]$ na 100 točk. Uporabite centralno diferenčno shemo drugega odvoda čez 5 točk, uteži izračunajte s pomočjo paketa findiff (funkcija findiff.coefficients).

Rezultat grafično prikažite (pazite na število robnih točk!).