Shared1. zadaća / 1.zadaca.ipynbOpen in CoCalc
Authors: Domagoj Bošnjak, Ivica Nakić
Views : 6
Description: 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 n×nn\times n linearnih jednadžbi Ax=bAx=b, potom zadaće izoliranja spektra proizvoljne matrice AA, 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. 2.2204e162.2204e^{-16} 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 10×1010\times 10 matricu i proivoljni vektor desne strane, b, došli do rješenja x.

In [ ]:
import numpy as np from numpy import linalg as lg A = np.random.rand(10,10) b = np.random.rand(10) if lg.det(A) != 0: x = np.dot(lg.inv(A),b) norma = lg.norm(b - np.dot(A,x)) print(norma)

Dakle, izračunati xx egzaktno rješava Ax=b+δbAx=b+\delta b, gdje je δbb21015\dfrac{||\delta b||}{||b||} \lessapprox 2\cdot 10^{-15} 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.

In [ ]:
import scipy.sparse.linalg as spla x, success = spla.gmres(A,b,tol=1e-10) print('Razlika izmedu egzaktnog rjesenja x i aprkosimacija za gmres metodu:{}, info:{}'.format(lg.norm(b - A@x), success)) x, success = spla.cg(A,b,tol=1e-10) print('Razlika izmedu egzaktnog rjesenja x i aprkosimacija za cg metodu:{}, info:{}'.format(lg.norm(b - A@x), success)) x, success = spla.bicg(A,b,tol=1e-10) print('Razlika izmedu egzaktnog rjesenja x i aprkosimacija za bicg metodu:{}, info:{}'.format(lg.norm(b - A@x), success))

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 A1A^{-1} tako da dodamo još jedan argument metodi cg(M=).

In [ ]:
x, success = spla.cg(A,b,tol=1e-10, M=lg.inv(A)) print('Razlika izmedu egzaktnog rjesenja x i aprkosimacija za cg metodu:{}, info:{}'.format(lg.norm(b - A@x), success))

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

In [ ]:
class gmres_counter(object): def __init__(self, disp=True): self._disp = disp self.niter = 0 def __call__(self, rk=None): self.niter += 1 if self._disp: print('iter %3i\trk = %s' % (self.niter, str(rk)))

Potom iskoristimo argument callback koji nam omogućuje da se nakon svake iteracije pozove proizvoljna funkcija.

In [ ]:
counter = gmres_counter() x, info = spla.gmres(A, b, callback=counter,tol=1e-10) print(counter.niter)

Numeričko rješavanje problema svojstvenih vrijednosti

Iz linerane algebre znamo da su nultočke karakterističnog polinoma svojstvene vrijednosti od linearnog operatora AA koji ima svoju matričnu reprezentaciju:

In [ ]:
print(A)

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

In [ ]:
from numpy.linalg import eig D,V = eig(A) print('D={}'.format(D)) print('V={}'.format(V)) print(lg.norm(np.dot(A[5,:],V[:,5]) - np.dot(D[5],V[:,5]))) print('Norma drugog stupca matrice D je {}'.format(lg.norm(V[:,2])))

Funckija dakle, vraća vektor svojstvenih vrijednosti DD i matricu svojstvenih (normaliziranih)vektora VV tako da je AV=VDA\cdot V = V\cdot D. Sada bismo htjeli doći do metoda koje bi računale istu stvar. Očito, za proizvoljnu matricu AA prvo moramo barem otprilike locirati nekako njen spektar. Jedan od načina su tzv. Geršgorinovi krugovi.

Teorem 1.\textbf{Teorem 1.} Neka je ACn×n\texttt{Neka je A$\in\mathbb{C}^{n\times n}$} i neka su λ1,,λn\lambda_1,\dots,\lambda_n njene svojstvene vrijednosti. Za i=1,,ni=1,\dots,n definirajmo Geršgorinove krugove: Gi={zC:zaiiρi},ρi=j=1jinaij\mathcal{G}_i = \{z\in\mathbb{C}:|z-a_{ii}|\leq \rho_i\},\quad \rho_i = \sum_{j=1\\j\neq i}^{n}|a_{ij}| Tada vrijedi:

\qquad \bulletSve svojstvene vrijednosti od AA su sadržane u uniji svih Geršgorinovih krugova, tj. σ(A)i=1nGi\sigma(A)\subseteq\bigcup_{i=1}^{n}\mathcal{G}_i

Nacrtajmo svojstvene vrijednosti naše matrice AA, tj. komponente vektora DD

In [ ]:
def draw_eigen(a): import matplotlib.pyplot as plt import numpy as np plt.grid(True) for x in range(len(a)): plt.plot([a[x].real],[a[x].imag], marker='o', markerfacecolor="yellow", markeredgewidth=4, markeredgecolor="blue", markersize=20) limit=np.max(np.ceil(np.absolute(a))) plt.xlim(-max(abs(D)) - 1,max(abs(D)) + 1) plt.ylim((-max(abs(D)) - 1,max(abs(D)) + 1)) plt.ylabel('Imaginary') plt.xlabel('Real') plt.show() draw_eigen(D)

Dakle, vidimo da su sve svojstvene vrijednosti negdje unutar kruga K(0,5)K(0,5)

Metoda potencija

Metoda potencija predstavlja najjednostavniju metodu za pronalaženje svojstvenih vrijednosti matrice AA. 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 A=uvA = uv^∗ , gdje su u,vCn, u,v0.u,v \in\mathcal{C}^n,\ u,v\neq0. Kako je Au=(vu)uAu = (v^∗u)u, vidimo da je λ1=vuλ_1 = v^∗u svojstvena vrijednost s pripadnim svojstvenim vektorom uu. Kako je ortogonalni komplement od L(v)\mathcal{L}(v) jezgra matrice AA, vidimo da je ostatak spektra nula, λ2==λn=0λ_2 = \dots = λ_n = 0.

Sada zamislimo da trebamo naći svojstveni vektor i svojstvenu vrijednost od AA i da znamo da je matrica AA ranga jedan, tj. da je oblika Au=(vu)uAu = (v^∗u)u, pri čemu su vektori u,vu,v nepoznati. Uzmimo proizvoljan vektor x0x\neq 0 i izračunajmo y=Ax=u(vx)y = Ax = u(v^∗x). Ako je y=0y = 0 onda je xx jedan svojstveni vektor svojstvene vrijednosti nula. Ako je y0y\neq 0, onda je yy kolinearan s vektorom uu (svojstvenim vektorom jedine netrivijalne svojstvene vrijednosti). Dakle je yy i sam svojstveni vektor, Ay=λ1yAy = λ_1y i svojstvena vrijednost λ1λ_1 se lako izračuna kao λ1=yAyyy.λ_1 = \dfrac{y^∗Ay}{y^∗ y}.

Napisimo funkciju koja za danu matricu A koja ima svojstvene vrijednosti λ1,λ2,,λnλ_1 , λ_2 ,\dots, λ_n , numerirane tako da je λ1>λ2λn|λ_1 | > |λ_2 | \geq \dots \geq |λ_n |, vrati aproksimaciju po modulu najvecu svojstvenu vrijednost.

In [ ]:
def pot_method(A, tol): VecRes = list() x0 = np.random.rand(10,1) ro = np.dot(x0.conj().T, np.dot(A,x0))/np.dot(x0.conj().T, x0) res = lg.norm(np.dot(A, x0) - ro*x0) VecRes.append(res) br = 1 while(res >= tol): br = br + 1; x = np.dot(A, x0)/lg.norm(np.dot(A,x0)) x0 = x ro = np.dot(x0.conj().T, np.dot(A,x0))/np.dot(x0.conj().T, x0) res = lg.norm(np.dot(A, x0) - ro*x0) VecRes.append(res) return [x, ro, VecRes, br, tol ] [x, ro, VecRes, br, tol ] = pot_method(A, 1e-10) print('Aproksimacija najvece po modulu svojstvene vrijednosti: {}'.format(ro[0,0])) print('Najveca po modulu svojstvena vrijednost: {}'.format(max(abs(D)))) print('Apsolutna razlika: {}'.format(lg.norm(ro - max(abs(D))))) print('Broj iteracija potreban za postizanje tolerancije({}): {}'.format(tol, br))

Prethodnu analizu smo napravili u specijalnom slučaju dijagonalizabile matrice AA u kojoj je dominantna po modulu svojstvena vrijednost strogo veća od ostalih, λ1>maxi=2:nλi|λ_1| > \max_{i=2:n}|λ_i|. Vidimo dakle da nam je potrebno svega 16 iteracija za postizanje točnosti reda 11. Nacrtajmo još vektor reziduala koji vraća metoda potencija.

In [ ]:
import matplotlib.pyplot as plt %matplotlib inline fig, ax = plt.subplots(figsize=(10, 4)) x = range(1, br+1) ax.plot(x, VecRes, 'r--o', markeredgewidth=2, markeredgecolor="blue") ax.set_yscale('log') ax.set_xlabel('x') ax.set_ylabel('VecRes') ax.set_title('Prikaz konvergencije reziduala u 0');

Dakle, za otkrivanje, tj. aproksimaciju neke svojstvene vrijednosti potreban nam je samo njen pripadni svojstveni vektor.

Definicija 2.\textbf{Definicija 2.} Neka je zadana matrica ACn×nA\in\mathbb{C}^{n\times n} i Cnx0\mathbb{C}^n\ni x\neq 0 definiramo Rayleighev kvocijent ρ=ρ(A,x)=xAxxx\rho = \rho(A,x) = \dfrac{x^*Ax}{x^*x}

Modificirajmo predloženu metodu da radi za bilo koju matricu dimenzije n×nn\times n

In [ ]:
from IPython.display import display from ipywidgets import interact, fixed import ipywidgets as widgets def pot_method_gen(n, tol): VecRes = list() A = np.random.rand(n,n) D,V = eig(A) x0 = np.random.rand(n,1) ro = np.dot(x0.conj().T, np.dot(A,x0))/np.dot(x0.conj().T, x0) res = lg.norm(np.dot(A, x0) - ro*x0) VecRes.append(res) br = 1 while(res >= tol): br = br + 1; x = np.dot(A, x0)/lg.norm(np.dot(A,x0)) x0 = x ro = np.dot(x0.conj().T, np.dot(A,x0))/np.dot(x0.conj().T, x0) res = lg.norm(np.dot(A, x0) - ro*x0) VecRes.append(res) fig, ax = plt.subplots(figsize=(10, 4)) domain = range(1, br+1) ax.plot(domain, VecRes, 'r--o', markeredgewidth=2, markeredgecolor="blue") ax.set_yscale('log') ax.set_xlabel('domain') ax.set_ylabel('VecRes') ax.set_title('Prikaz konvergencije reziduala u 0'); ax.grid('on') return interact(pot_method_gen,n=(1,10), tol=fixed(1e-15));

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:

Objašnjenje involviranih objekata u nasoj zadaći:

f: I×RdRdf:\ \mathcal{I}\times\mathbb{R}^d\longrightarrow\mathbb{R}^d gdje je IR\mathcal{I}\subseteq\mathbb{R} otvoreni interval, t0I,d1t_0\in\mathcal{I},d\geq 1 i y0Rdy_0\in\mathbb{R}^d te očito onda je y:IRdy:\mathcal{I}\longrightarrow\mathbb{R}^d. Ako je d2d\geq 2 onda govorimo o sustavu običnih diferencijalnih jednadžbi i to možemo prikazati na sljedeći način(po komponentama):

(y1(t)y2(t)yd(t))=(f1(t,y1(t),y2(t),,yd(t))f2(t,y1(t),y2(t),,yd(t))fd(t,y1(t),y2(t),,yd(t))),(y1(t0)y2(t0)yd(t0))=((y0)1(y0)2(y0)d) \begin{pmatrix} y_1'(t)\\ y_2'(t)\\ \vdots\\ y_d'(t)\\ \end{pmatrix} = \begin{pmatrix} f_1(t,y_1(t),y_2(t),\dots,y_d(t))\\ f_2(t,y_1(t),y_2(t),\dots,y_d(t))\\ \vdots\\ f_d(t,y_1(t),y_2(t),\dots,y_d(t))\\ \end{pmatrix} ,\quad \begin{pmatrix} y_1(t_0)\\ y_2(t_0)\\ \vdots\\ y_d(t_0)\\ \end{pmatrix} = \begin{pmatrix} (y_0)_{1}\\ (y_0)_{2}\\ \vdots\\ (y_0)_{d}\\ \end{pmatrix} I to zovemo inicijalni problem za sustav običnih diferencijalnih jednadžbi.Ako je

f(t,x1,,xd)=A(t)(x1x2xd)+(b1(t)b2(t)bd(t))f(t,x_1,\dots,x_d) = A(t) \begin{pmatrix} x_1\\ x_2\\ \vdots\\ x_d \end{pmatrix} + \begin{pmatrix} b_1(t)\\ b_2(t)\\ \vdots\\ b_d(t) \end{pmatrix} Onda kažemo da je sustav ODJ linearan.

Eulerova metoda

Eulerova metoda je najjednostavniji primjer jednokoracˇne eksplicitne metode;\texttt{jednokoračne eksplicitne metode;} vrijednost yi+1y_{i+1} je dobivena eksplicitnim izrazom koji koristi informaciju samo iz koraka tit_i. 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: y(0)=1y(0) = 1 i y(10)=2y(10) = 2

Egzatkno rješenje je dano s:

In [ ]:
import math as mt def funname(t, y): f = [y[0] - 2*y[1] + 4*mt.cos(t) - 2*mt.sin(t), 3*y[0] - 4*y[1] + 5*mt.cos(t) - 5*mt.sin(t)] return f def euler_method(funname, y0, a, b, n): h = (b-a)/n y = np.random.rand(n+1,2) y[0,:] = y0 egz = np.random.rand(n+1,2) egz[0,:] = y0 t = np.random.rand(n+1,1) t[0] = a; for i in range(0,n): t[i+1] = t[i] + h egz[i+1, 0] = mt.cos(t[i+1]) + mt.sin(t[i+1]) egz[i+1, 1] = 2*mt.cos(t[i+1]) rj = list() rj.append(h*funname(t[i], y[i,:])[0]) rj.append(h*funname(t[i], y[i,:])[1]) y[i+1,:] = y[i,:] + rj error = list() err = max(abs(y[:,0] - egz[:,0])) error.append(err) err = max(abs(y[:,1] - egz[:,1])) error.append(err) fig, ax = plt.subplots(figsize=(20, 6)) ax.plot(t, egz[:,0], '-b') ax.plot(t, egz[:,1], '-g') ax.plot(t, y[:,0], '--r') ax.plot(t, y[:,1], '--k') ax.set_xlabel('t') ax.legend(["egzaktno rjesenje y_1", "egzaktno rjesenje y_2", "aproksimacija rjesenja y_1", "aproksimacija rjesenja y_2"]); return [y,t, error] [y,t, error] = euler_method(funname, [1,2], 0, 10, 100) print('Globalna greska diskretizacije Eulerove metode za y(1) je :{}'.format(error[0])) print('Globalna greska diskretizacije Eulerove metode za y(2) je :{}'.format(error[1]))

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\cdot0.01(C>0)(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

In [ ]:
from scipy.integrate import odeint as ont def ode_sys(y,t): y1,y2 = y dydt = [y1 - 2*y2 + 4*mt.cos(t)-2*mt.sin(t), 3*y1 - 4*y2 + 5*mt.cos(t) - 5*mt.sin(t)] return dydt y0 = [1,2] t = np.linspace(0,10,100) sol = ont(ode_sys, y0, t) import matplotlib.pyplot as plt plt.plot(t, sol[:, 0], 'b', label='y1(t)') plt.plot(t, sol[:, 1], 'g', label='y2(t)') plt.legend(loc='best') plt.xlabel('t') plt.grid() plt.show()

Implementirajmo sada Runge-Kutta metoda sa stupnjem i redom jednak 4 za istu zadacu ()(\bigstar). Ona je i dalje predstavnik eksplicitne metode, ali sada se traženi integral aproksimira tako da se između koraka diskretizacije ti,ti+1t_i,t_{i+1} 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

In [ ]:
def RK_4_classic(funname, y0, a, b, n): h = (b-a)/n y = np.random.rand(n+1,2) y[0,:] = y0 egz = np.random.rand(n+1,2) egz[0,:] = y0 t = np.random.rand(n+1,1) t[0] = a; for i in range(0,n): t[i+1] = t[i] + h egz[i+1, 0] = mt.cos(t[i+1]) + mt.sin(t[i+1]) egz[i+1, 1] = 2*mt.cos(t[i+1]) k1 = funname(t[i],y[i,:]) k2 = funname(t[i] + (1/2)*h, y[i,:] + [h*(1/2)*z for z in k1]) k3 = funname(t[i] + (1/2)*h, y[i,:] + [h*(1/2)*z for z in k2]) k4 = funname(t[i] + h, y[i,:] + [h*z for z in k3]) k2 = [2*x for x in k2] k2 = [3*x for x in k3] y[i+1,:] = y[i,:] + [(1/6)*h*z for z in [a1 + a2 + a3 + a4 for a1,a2,a3,a4 in zip(k1,k2,k3,k4)]] error = list() err = max(abs(y[:,0] - egz[:,0])) error.append(err) err = max(abs(y[:,1] - egz[:,1])) error.append(err) fig, ax = plt.subplots(figsize=(20, 6)) ax.plot(t, egz[:,0], '-b') ax.plot(t, egz[:,1], '-g') ax.plot(t, y[:,0], '--r') ax.plot(t, y[:,1], '--k') ax.set_xlabel('t') ax.legend(["egzaktno rjesenje y_1", "egzaktno rjesenje y_2", "aproksimacija rjesenja y_1", "aproksimacija rjesenja y_2"]); return [y,t, error] [y,t, error] = RK_4_classic(funname, [1,2], 0, 10, 100) print('Globalna greska diskretizacije Runge-Kutta metode za y(1) je :{}'.format(error[0])) print('Globalna greska diskretizacije Runge-Kutta za y(2) je :{}'.format(error[1]))

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

In [ ]:
def RK_4_classic(funname, y0, a, b, n): h = (b-a)/n y = np.random.rand(n+1,2) y[0,:] = y0 egz = np.random.rand(n+1,2) egz[0,:] = y0 t = np.random.rand(n+1,1) t[0] = a; for i in range(0,n): t[i+1] = t[i] + h egz[i+1, 0] = mt.cos(t[i+1]) + mt.sin(t[i+1]) egz[i+1, 1] = 2*mt.cos(t[i+1]) k1 = funname(t[i],y[i,:]) k2 = funname(t[i] + (1/2)*h, y[i,:] + [h*(1/2)*z for z in k1]) k3 = funname(t[i] + (1/2)*h, y[i,:] + [h*(1/2)*z for z in k2]) k4 = funname(t[i] + h, y[i,:] + [h*z for z in k3]) k2 = [2*x for x in k2] k2 = [3*x for x in k3] y[i+1,:] = y[i,:] + [(1/6)*h*z for z in [a1 + a2 + a3 + a4 for a1,a2,a3,a4 in zip(k1,k2,k3,k4)]] error = list() err = max(abs(y[:,0] - egz[:,0])) error.append(err) err = max(abs(y[:,1] - egz[:,1])) error.append(err) fig, ax = plt.subplots(figsize=(20, 6)) ax.plot(t, egz[:,0], '-b') ax.plot(t, egz[:,1], '-g') ax.plot(t, y[:,0], '--r') ax.plot(t, y[:,1], '--k') ax.set_xlabel('t') ax.legend(["egzaktno rjesenje y_1", "egzaktno rjesenje y_2", "aproksimacija rjesenja y_1", "aproksimacija rjesenja y_2"]); print('Globalna greska diskretizacije Runge-Kutta metode za y(1) je :{}'.format(error[0])) print('Globalna greska diskretizacije Runge-Kutta za y(2) je :{}'.format(error[1])) return interact(RK_4_classic,n=(1,10000), funname=fixed(funname), y0=fixed([1,2]), a=fixed(0), b=fixed(10));

Rijesimo opet našu zadaću ()(\bigstar), sada koristeći klasu ode iz modula scipy.integrate

In [ ]:
import math as mt from scipy.integrate import ode import matplotlib.pyplot as plt def f(t, y): return [y[0] - 2*y[1] + 4*mt.cos(t) - 2*mt.sin(t), 3*y[0] - 4*y[1] + 5*mt.cos(t) - 5*mt.sin(t)] y0 = [1,2] t0 = 0 t1= 10 N = 1000 y = np.empty((N+1,2)) egz = np.empty((N+1,2)) y[0] = y0 egz[0] = y0 k = 1 t = np.linspace(t0, t1, N+1) r = ode(f).set_integrator('dopri5') r.set_initial_value(y0, t0) while r.successful() and r.t < t1: r.integrate(t[k]) y[k] = r.y k += 1 for i in range(0,N): egz[i+1, 0] = mt.cos(t[i+1]) + mt.sin(t[i+1]) egz[i+1, 1] = 2*mt.cos(t[i+1]) error = list() err = max(abs(y[:,0] - egz[:,0])) error.append(err) err = max(abs(y[:,1] - egz[:,1])) error.append(err) fig, ax = plt.subplots(figsize=(20, 6)) ax.plot(t, egz[:,0], '-b') ax.plot(t, egz[:,1], '-g') ax.plot(t, y[:,0], '--r') ax.plot(t, y[:,1], '--k') ax.set_xlabel('t') ax.legend(["egzaktno rjesenje y_1", "egzaktno rjesenje y_2", "aproksimacija rjesenja y_1", "aproksimacija rjesenja y_2"]); print('Globalna greska diskretizacije Runge-Kutta metode za y(1) je :{}'.format(error[0])) print('Globalna greska diskretizacije Runge-Kutta za y(2) je :{}'.format(error[1]))

Analiza podataka

%Nisam uspio instalirati pakete pa sam ostao na ovom jednostavnom primjeru

In [ ]:
import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.basemap import Basemap %matplotlib inline themap = Basemap(projection='gall') themap.drawcoastlines() themap.drawcountries() themap.drawstates() themap.fillcontinents(color = 'gainsboro') themap.drawmapboundary(fill_color='blue');

Demografija

Primjer nekoliko metoda unutar modula pandas. Podaci su preuzeti sa Link.

In [ ]:
import pandas as pd names = ['Area', 'Year','Area', 'Sex', 'Age','School','attendance','Record','Type','Reliability','Source','Year','Value'] data1 = pd.read_csv('albania5_24.csv', delimiter=',', names=names) data2 = pd.read_csv('population_sex.csv', delimiter=',', names=names) data3 = pd.read_csv('albania_15_and_over.csv', delimiter=',', names=names) print(data1.shape) print(data2.shape) print(data3.shape) print(data1.head(5)) print(data1.first_valid_index()) print(data2.get_values()) print(data3.keys()) data3_dict = data3.to_dict() print(data2.info())
In [ ]: