1.zadaca.ipynb

Kernel: Python 3 (old Anaconda 3)

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\times n$ linearnih jednadžbi $Ax=b$, potom zadaće izoliranja spektra proizvoljne matrice $A$, 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.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\times 10$ matricu i proivoljni vektor desne strane, b, došli do rješenja x.

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 $x$ egzaktno rješava $Ax=b+\delta b$, gdje je $\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.

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

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

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.

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 $A$ koji ima svoju matričnu reprezentaciju:

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

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 $D$ i matricu svojstvenih (normaliziranih)vektora $V$ tako da je $A\cdot V = V\cdot D$. Sada bismo htjeli doći do metoda koje bi računale istu stvar. Očito, za proizvoljnu matricu $A$ prvo moramo barem otprilike locirati nekako njen spektar. Jedan od načina su tzv. Geršgorinovi krugovi.

$\textbf{Teorem 1.}$ 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 $\lambda_1,\dots,\lambda_n$ njene svojstvene vrijednosti. Za $i=1,\dots,n$ definirajmo Geršgorinove krugove: $$\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$ $\bullet$Sve svojstvene vrijednosti od $A$ su sadržane u uniji svih Geršgorinovih krugova, tj. $\sigma(A)\subseteq\bigcup_{i=1}^{n}\mathcal{G}_i$

Nacrtajmo svojstvene vrijednosti naše matrice $A$, tj. komponente vektora $D$

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)$

Metoda potencija

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

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

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

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 $A$ u kojoj je dominantna po modulu svojstvena vrijednost strogo veća od ostalih, $|λ_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.

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.

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

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

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: ParseError: KaTeX parse error: {align*} can be used only in display mode.

Objašnjenje involviranih objekata u nasoj zadaći:

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

$$\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,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 $\texttt{jednokoračne eksplicitne metode;}$ vrijednost $y_{i+1}$ je dobivena eksplicitnim izrazom koji koristi informaciju samo iz koraka $t_i$. Eulerova metoda predstavlja paradigmu aproksimacije traženoga integrala pomoću pravokutnika

Definirajmo funkciju euler_method koja će riješiti sljedeću zadaću:

$$\begin{align*} (\bigstar)=\begin{cases} y(1)' &= y(1) - 2y(2) + 4\cos(t) - 2\sin(t)\\ y(2)' &= 3y(1) - 4y(2) + 5\cos(t) - 5\sin(t) \end{cases} \end{align*}$$

Uz početne uvjete: $y(0) = 1$ i $y(10) = 2$

Egzatkno rješenje je dano s: $$\begin{align*} y(1) &= \cos(t) + \sin(t)\\ y(2) &= 2\cos(t) \end{align*}$$

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$\cdot$0.01**$(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

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 $t_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

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 $y_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.

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

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

$$

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.

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())