Shared1. zadaća / 1.zadaca.ipynbOpen in CoCalc
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.

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)
1.44068549693e-15

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.

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))
Razlika izmedu egzaktnog rjesenja x i aprkosimacija za gmres metodu:1.121270291988505e-15, info:0 Razlika izmedu egzaktnog rjesenja x i aprkosimacija za cg metodu:94537.85488755246, info:100 Razlika izmedu egzaktnog rjesenja x i aprkosimacija za bicg metodu:5.351174074293835e-12, info:0

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=).

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))
Razlika izmedu egzaktnog rjesenja x i aprkosimacija za cg metodu:1.8975605279926845e-14, info:0

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)
iter 1 rk = 0.5674607588331203 iter 2 rk = 0.5665774550546372 iter 3 rk = 0.5623937145628866 iter 4 rk = 0.4645746507839533 iter 5 rk = 0.4624418315395474 iter 6 rk = 0.4615021046395359 iter 7 rk = 0.4590025930410239 iter 8 rk = 0.1556972760824254 iter 9 rk = 0.14423589605795206 iter 10 rk = 2.5354718799162735e-16 10

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:

print(A)
[[ 0.02419054 0.50955696 0.43238954 0.1551312 0.24033346 0.20985601 0.80275121 0.12786896 0.86377098 0.29284547] [ 0.98581118 0.76897724 0.6308045 0.14379682 0.94893868 0.04550622 0.05553063 0.80723444 0.99694585 0.16138347] [ 0.1010132 0.6571242 0.55580273 0.27013769 0.62739755 0.18903114 0.21609265 0.77806164 0.85823304 0.30802321] [ 0.49605583 0.38034957 0.77293409 0.08284963 0.91113489 0.63809338 0.14745068 0.33166553 0.60522122 0.31491286] [ 0.5497267 0.22769256 0.95767161 0.93271692 0.86254442 0.99026325 0.74406014 0.43240183 0.63955887 0.968017 ] [ 0.67969345 0.71005073 0.84241557 0.20623018 0.37569083 0.29954106 0.20735126 0.46521642 0.46082809 0.91814322] [ 0.70439017 0.85534618 0.58196472 0.10551355 0.22452321 0.8840903 0.43190607 0.19117882 0.28224431 0.82637858] [ 0.74359067 0.03934782 0.64830243 0.42989314 0.59717646 0.32510503 0.67254882 0.56835037 0.10099016 0.03825684] [ 0.84969367 0.03731904 0.57043543 0.71846879 0.56839042 0.31522219 0.28181657 0.01008406 0.26005889 0.42355342] [ 0.04762778 0.95931947 0.3995724 0.82914664 0.53916876 0.34707601 0.6287983 0.07033233 0.37015683 0.39594176]]

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])))
D=[ 4.92148928+0.j 0.16623371+1.01396985j 0.16623371-1.01396985j 0.59038501+0.16892665j 0.59038501-0.16892665j -0.89442992+0.j -0.51642751+0.46596629j -0.51642751-0.46596629j -0.35930913+0.j 0.10203006+0.j ] V=[[ 2.30105891e-01+0.j 1.20466068e-01-0.03288589j 1.20466068e-01+0.03288589j -7.68985830e-02+0.31668283j -7.68985830e-02-0.31668283j 4.77517675e-01+0.j 5.21379726e-02-0.00208301j 5.21379726e-02+0.00208301j -2.29513302e-02+0.j 2.14453295e-02+0.j ] [ 3.47643763e-01+0.j -5.31914843e-01+0.j -5.31914843e-01-0.j 4.39029536e-02+0.09382129j 4.39029536e-02-0.09382129j -3.19349346e-02+0.j -2.85460628e-01+0.0626842j -2.85460628e-01-0.0626842j -6.27958368e-03+0.j -4.98706686e-02+0.j ] [ 2.91389861e-01+0.j -3.15571643e-01+0.0977758j -3.15571643e-01-0.0977758j 1.61874924e-01-0.25346093j 1.61874924e-01+0.25346093j 4.56323232e-01+0.j -1.71961361e-01-0.02715583j -1.71961361e-01+0.02715583j -5.00644002e-01+0.j -5.76002226e-01+0.j ] [ 3.06736856e-01+0.j -7.18960581e-02-0.05544169j -7.18960581e-02+0.05544169j -1.18620027e-01-0.19400459j -1.18620027e-01+0.19400459j 9.48023062e-03+0.j -1.68946891e-01+0.39083136j -1.68946891e-01-0.39083136j 3.29010342e-01+0.j -2.41365287e-01+0.j ] [ 4.66871992e-01+0.j 2.96007930e-01-0.23506717j 2.96007930e-01+0.23506717j -1.40041315e-01-0.20899617j -1.40041315e-01+0.20899617j 4.20196129e-02+0.j 2.57136549e-04-0.12035987j 2.57136549e-04+0.12035987j -3.44827592e-01+0.j 5.15431740e-01+0.j ] [ 3.18828064e-01+0.j 5.03751524e-02+0.25945411j 5.03751524e-02-0.25945411j 7.95439217e-02+0.00640567j 7.95439217e-02-0.00640567j -3.92299280e-01+0.j -3.42278789e-01-0.23224034j -3.42278789e-01+0.23224034j 4.78914209e-01+0.j -1.64240929e-01+0.j ] [ 3.17045498e-01+0.j 3.57233299e-01+0.25644268j 3.57233299e-01-0.25644268j 1.72923802e-01+0.33933925j 1.72923802e-01-0.33933925j -1.40062407e-01+0.j -5.43726121e-02+0.1771419j -5.43726121e-02-0.1771419j -1.79858641e-01+0.j 3.00946806e-01+0.j ] [ 2.61562113e-01+0.j 1.03716167e-02-0.27877793j 1.03716167e-02+0.27877793j 5.97191523e-01+0.j 5.97191523e-01-0.j -2.80958251e-01+0.j 2.43726745e-01+0.00421205j 2.43726745e-01-0.00421205j 3.68690194e-01+0.j 1.47477182e-02+0.j ] [ 2.53778987e-01+0.j 8.66509154e-02-0.17676513j 8.66509154e-02+0.17676513j -3.68327805e-01+0.05054872j -3.68327805e-01-0.05054872j -5.25073772e-01+0.j 1.41457734e-01-0.13948685j 1.41457734e-01+0.13948685j 3.34688683e-01+0.j -1.71708505e-01+0.j ] [ 3.06989236e-01+0.j 8.34789654e-02+0.23329611j 8.34789654e-02-0.23329611j -1.74359773e-01-0.03068659j -1.74359773e-01+0.03068659j 1.80871245e-01+0.j 6.17624464e-01+0.j 6.17624464e-01-0.j -1.09446743e-01+0.j 4.40570086e-01+0.j ]] 1.37952226836 Norma drugog stupca matrice D je 1.0

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.} i neka su λ1,,λn\lambda_1,\dots,\lambda_n njene svojstvene vrijednosti. Za i=1,,ni=1,\dots,n definirajmo Geršgorinove krugove: 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

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 , gdje su u,vCn, u,v0.u,v \in\mathcal{C}^n,\ u,v\neq0. Kako je , vidimo da je 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 , pri čemu su vektori u,vu,v nepoznati. Uzmimo proizvoljan vektor x0x\neq 0 i izračunajmo . 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

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.

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))
Aproksimacija najvece po modulu svojstvene vrijednosti: 4.9214892836857045 Najveca po modulu svojstvena vrijednost: 4.921489283686206 Apsolutna razlika: 5.018208071305708e-13 Broj iteracija potreban za postizanje tolerancije(1e-10): 16

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.

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

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));
Unsupported message: {"data":{"application/vnd.jupyter.widget-view+json":{"model_id":"d187bf15041848249ba32e74f9489e7d"}},"execution_count":29,"metadata":{},"output_type":"execute_result"}

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 jednokorane eksplicitne metode;cˇ\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:

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]))
Globalna greska diskretizacije Eulerove metode za y(1) je :0.05889826324201239 Globalna greska diskretizacije Eulerove metode za y(2) je :0.036670785225934255

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

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

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]))
Globalna greska diskretizacije Runge-Kutta metode za y(1) je :0.0008102189436380858 Globalna greska diskretizacije Runge-Kutta za y(2) je :0.0016793178043463808

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.

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));
Unsupported message: {"data":{"application/vnd.jupyter.widget-view+json":{"model_id":"e9b5dfd9498b4cb2a25dac260e461dcf"}},"execution_count":33,"metadata":{},"output_type":"execute_result"}

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]))
Globalna greska diskretizacije Runge-Kutta metode za y(1) je :5.146993942162226e-13 Globalna greska diskretizacije Runge-Kutta za y(2) je :8.038014698286133e-13

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())
/projects/anaconda3/lib/python3.5/site-packages/IPython/core/interactiveshell.py:2683: DtypeWarning: Columns (8,9) have mixed types. Specify dtype option on import or set low_memory=False. interactivity=interactivity, compiler=compiler, result=result)
(86133, 13) (78248, 13) (98060, 13) Area Year Area.1 Sex Age School \ 0 Country or Area Year Area Sex Age School attendance 1 Albania 2011 Total Both Sexes 6 Total 2 Albania 2011 Total Both Sexes 6 Attending school 3 Albania 2011 Total Both Sexes 6 Not attending school 4 Albania 2011 Total Both Sexes 6 - 23 Total attendance Record \ 0 Record Type Reliability 1 Census - de jure - complete tabulation Final figure, complete 2 Census - de jure - complete tabulation Final figure, complete 3 Census - de jure - complete tabulation Final figure, complete 4 Census - de jure - complete tabulation Final figure, complete Type Reliability Source Year.1 Value 0 Source Year Value Value Footnotes NaN NaN 1 2013 33833 NaN NaN NaN 2 2013 33833 NaN NaN NaN 3 2013 0 NaN NaN NaN 4 2013 827512 NaN NaN NaN 0 [['Country or Area' 'Year' 'Area' ..., 'Value Footnotes' nan nan] ['Albania' '2011' 'Total' ..., nan nan nan] ['Albania' '2011' 'Total' ..., nan nan nan] ..., ['49' 'Excluding population enumerated in hotels.' nan ..., nan nan nan] ['50' 'Excluding homeless persons.' nan ..., nan nan nan] ['51' 'Excluding Indian jungle population.' nan ..., nan nan nan]] Index(['Area', 'Year', 'Area.1', 'Sex', 'Age', 'School', 'attendance', 'Record', 'Type', 'Reliability', 'Source', 'Year.1', 'Value'], dtype='object') <class 'pandas.core.frame.DataFrame'> RangeIndex: 78248 entries, 0 to 78247 Data columns (total 13 columns): Area 78248 non-null object Year 78248 non-null object Area.1 78196 non-null object Sex 78196 non-null object Age 78196 non-null object School 78196 non-null object attendance 78196 non-null object Record 78196 non-null object Type 78196 non-null object Reliability 78196 non-null object Source 16569 non-null object Year.1 0 non-null float64 Value 0 non-null float64 dtypes: float64(2), object(11) memory usage: 7.8+ MB None