Jupyter notebook 1. zadaća/1.zadaca.ipynb
Numerička analiza
Numerička analiza je dio matematike koja proučava metode kojima se pokušava aproksimirati rješenje dane zadaće. Involvirane zadaće se najčešće odnose na rješavanje sustava linearnih jednadžbi , potom zadaće izoliranja spektra proizvoljne matrice , numeričke aproksimacije sustava ODJ-a itd. Naglasak je na aproksimaciji jer do egzaktnih rješenja u praksi dolazimo vrlo rijetko. Jedan od razloga je pogreška strojne aritmetike, tj. neminovno je da će rčunalo čak pri zbrajanju dvaju brojeva pogriješiti, tj. reprezentirano rješenje neće biti egzatkno, ali od njega će biti udaljeno "jako malo", ovisnosti o preciznosti samoga stroja, npr. daleko. Očekivano, iako su teorijske ideje rješavanja nekih zadaća bile poznate i prije 18. stoljeća, samo područje numeričke analize je doživjelo pravi procvat razvojem računala, dakle drugom polovicom 20. stoljeca i danas predstavlja vrlo bitan kotačić u raznim granama industrije, građevini, hidrometeorološkim zavodima, medicini itd.
Numeričke metode za rješavanje sustava Ax=b
Iz linearne algebre znamo da zadani sustav ima rješenje ako i samo ako je determinanta matrice A različita od 0, tj. akko je matrica regularna. Provjerimo kako bismo za proizvoljnu matricu i proivoljni vektor desne strane, b, došli do rješenja x.
Dakle, izračunati egzaktno rješava , gdje je s čim bi više manje svatko bio zadovoljan. Izlažem metodu gmres, cg, bicg kojom bi pomaknuli rješavanje sustava pomoću operatora \ iz matlaba ili implementirane metode inv kroz određene matematičke alate, npr. Krilovljeve potprostore.
Vidimo da su metode gmres i bicg uredno iskonvergirale (jer je info==0), i da je aproksimacija vrlo blizu egazktnom rješenju. S druge strane, očitavamo da je ukupan broj iteracija metode cg jednak 100, a iz normi reziduala očito metoda nije iskonvergirala. Popravimo ju uz "too good to be true" prekondicioner tako da dodamo još jedan argument metodi cg(M=).
Očito ovako jaki kondicioner je pomogao funkciji da postigne konvergenciju. Također se iz dokumentacije Link1, Link2, Link3 može vidjeti da metodama možemo poslati maksimalan broj iteracija u slučaju da se ne postigne konvergencija, ali nigdje nema podatka o ukupnom broju iteracija dok metoda ne iskonvergira. To možemo dodati na sljedeći nacin: prvo definirajmo klasu gmres_counter
Potom iskoristimo argument callback koji nam omogućuje da se nakon svake iteracije pozove proizvoljna funkcija.
Numeričko rješavanje problema svojstvenih vrijednosti
Iz linerane algebre znamo da su nultočke karakterističnog polinoma svojstvene vrijednosti od linearnog operatora koji ima svoju matričnu reprezentaciju:
Općenito, izračun spektra prozivoljne matrice je skup i kroz vrijeme su razvijene razne numeričke metode koje su generalno manje ili više primijenjive na gotovo sve matrice. Za početak iskoristimo implementrianu funkciju eig i norm iz paketa numpy
Funckija dakle, vraća vektor svojstvenih vrijednosti i matricu svojstvenih (normaliziranih)vektora tako da je . Sada bismo htjeli doći do metoda koje bi računale istu stvar. Očito, za proizvoljnu matricu prvo moramo barem otprilike locirati nekako njen spektar. Jedan od načina su tzv. Geršgorinovi krugovi.
ParseError: KaTeX parse error: Unexpected end of input in a macro argument, expected '}' at end of input: …exttt{Neka je A\in\mathbb{C}^{n\times n}ParseError: KaTeX parse error: Expected 'EOF', got '}' at position 1: }̲ i neka su njene svojstvene vrijednosti. Za definirajmo Geršgorinove krugove: Tada vrijedi:
Sve svojstvene vrijednosti od su sadržane u uniji svih Geršgorinovih krugova, tj.
Nacrtajmo svojstvene vrijednosti naše matrice , tj. komponente vektora
Dakle, vidimo da su sve svojstvene vrijednosti negdje unutar kruga
Metoda potencija
Metoda potencija predstavlja najjednostavniju metodu za pronalaženje svojstvenih vrijednosti matrice . Jednostavnom primjenom matrice A na polazni vektor x se generira niz vektora koji pod određenim uvjetima daje informaciju o najvećoj po modulu svojstvenoj vrijednosti i pripadnom vektoru.Nekoliko je elementarnih ideja koje motiviraju metodu potencija. Za početak, promotrimo problem svojstvenih vrijednosti za matricu , gdje su Kako je , vidimo da je svojstvena vrijednost s pripadnim svojstvenim vektorom . Kako je ortogonalni komplement od jezgra matrice , vidimo da je ostatak spektra nula, .
Sada zamislimo da trebamo naći svojstveni vektor i svojstvenu vrijednost od i da znamo da je matrica ranga jedan, tj. da je oblika , pri čemu su vektori nepoznati. Uzmimo proizvoljan vektor i izračunajmo . Ako je onda je jedan svojstveni vektor svojstvene vrijednosti nula. Ako je , onda je kolinearan s vektorom (svojstvenim vektorom jedine netrivijalne svojstvene vrijednosti). Dakle je i sam svojstveni vektor, i svojstvena vrijednost se lako izračuna kao
Napisimo funkciju koja za danu matricu A koja ima svojstvene vrijednosti , numerirane tako da je , vrati aproksimaciju po modulu najvecu svojstvenu vrijednost.
Prethodnu analizu smo napravili u specijalnom slučaju dijagonalizabile matrice u kojoj je dominantna po modulu svojstvena vrijednost strogo veća od ostalih, . Vidimo dakle da nam je potrebno svega 16 iteracija za postizanje točnosti reda 11. Nacrtajmo još vektor reziduala koji vraća metoda potencija.
Dakle, za otkrivanje, tj. aproksimaciju neke svojstvene vrijednosti potreban nam je samo njen pripadni svojstveni vektor.
Neka je zadana matrica i definiramo Rayleighev kvocijent
Modificirajmo predloženu metodu da radi za bilo koju matricu dimenzije
Numeričko rješavanje inicijalnog problema za ODJ
U nastavku se okrećemo, mozda čak i "najkorisnijem" području primjene numeričke analize, a to je rješavanje (sustava) ODJ-a i PDJ-a. Pokazat ću neke metode za aproksimaciju inicijalnoga problema za sustav ODJ-a
Zadaća koja je pred nama ima sljedeći oblik: ParseError: KaTeX parse error: {align*} can be used only in display mode.
Objašnjenje involviranih objekata u nasoj zadaći:
gdje je otvoreni interval, i te očito onda je . Ako je onda govorimo o sustavu običnih diferencijalnih jednadžbi i to možemo prikazati na sljedeći način(po komponentama):
I to zovemo inicijalni problem za sustav običnih diferencijalnih jednadžbi.Ako je
Onda kažemo da je sustav ODJ linearan.
Eulerova metoda
Eulerova metoda je najjednostavniji primjer vrijednost je dobivena eksplicitnim izrazom koji koristi informaciju samo iz koraka . Eulerova metoda predstavlja paradigmu aproksimacije traženoga integrala pomoću pravokutnika
Definirajmo funkciju euler_method koja će riješiti sljedeću zadaću:
Uz početne uvjete: i
Egzatkno rješenje je dano s:
Sa prikazanoga grafa vidimo da nam se rješenja poklapaju naizgled skoro savršeno, ali to nije tako. Zapravo pošto je Eulerova metoda primjer najjednostavnije metode za aproksimaciju ODJ-a ona ima najveću lokalnu pogrešku diskretizacije; za h=0.1, lok. pogreška će biti manja od **C$\cdot(C>0)$, dakle Eulerova metoda je reda 1. Naravno, postoje sofisticiranije metode od Eulerove koje imaju veći red (2,3,4,5,..) kao što su implicitna trapezna, eksplicitna Runge Kutta metoda koju ću poslije posebno obraditi, eksplicitne/implicitne višekoračne metode itd.
Riješimo za početak, sada našu zadacu pomoću implementirane funkcije odeint iz modula scipy.integrate
Implementirajmo sada Runge-Kutta metoda sa stupnjem i redom jednak 4 za istu zadacu . Ona je i dalje predstavnik eksplicitne metode, ali sada se traženi integral aproksimira tako da se između koraka diskretizacije ubaci po nekoliko čvorova, u našem slučaju 4, i proizvoljni parametri se uzmu tako da je lokalna pogreška diskretizacije minimalna. Tražene parametre sam uzeo sa sljedeće stranice : wiki
Vidimo sa usporedbi grafova da je Runga-Kutte metoda "jača" nego Eulerova, a iz prikazanih globalnih grešaka diskretizacije vidimo i koliko. Mane ove metode su što metoda generalno može biti skupa, u značenju da za stupanj=6 morat ćemo za jedan korak aproksimacije 6 puta izvrijednjavati funkciju desne strane; također, kao i svaka druga eksplicitna metoda niti Runge-Kutta se ne može koristiti za riješavanje tzv. "krutih sistema ODJ". Unatoč, metoda se vrlo često koristi, pa je i implementirana u modulu scipy što izlažem u sljedećoj ćeliji, ali za početak, postavimo da korisnik može sam birati broj točaka diskretizacije i pomoću istoga možemo analizirati kako to utječe na globalne greške diskretizacije.
Rijesimo opet našu zadaću , sada koristeći klasu ode iz modula scipy.integrate
Analiza podataka
Demografija
Primjer nekoliko metoda unutar modula pandas. Podaci su preuzeti sa Link.