Numerička integracija¶

Prvo implementiramo pravokutnu formulu

In [ ]:
import numpy as np

def pravokutna(f,a,b,N=50):
    h =  (b-a)/N
    x = np.linspace(a+h/2,b-h/2,N) 
    y = f(x)
    h = (b - a)/N
    Mn = np.sum(h*y)
    return Mn

Ovo je primjer implementacije trapezne formule

In [3]:
def trapezna(f,a,b,N=50):
    x = np.linspace(a,b,N+1) 
    y = f(x)
    y_right = y[0:N] # lijevi rub intervala
    y_left = y[1:N+1] # desni rub intervala
    h = (b - a)/N
    Tn = np.sum(h*(y_right + y_left)/2)
    return Tn

Ovo je primjer implementacije Simpsonove formule

In [4]:
def Simpson(f,a,b,N=50):
    h =  (b-a)/N
    x = np.linspace(a,b,N+1) 
    y = f(x)
    y[1:N] = 2*y[1:N]
    Sn = h*np.sum(y)/6
    x = np.linspace(a+h/2,b-h/2,N) # desni rub intervala
    y=f(x)
    Sn= Sn+2*h*np.sum(y)/3
    return Sn

Primijenit ćemo navedene formule na funkciju $f(x)=xe^{-x}\cos(2x)$ na segmentu $[0,2\pi]$

In [13]:
import numpy as np
from numpy import exp
from numpy import cos
from numpy import pi
print("Pravokutna formula daje:" , pravokutna(lambda x : x*exp(-x)*cos(2*x),0,2*pi,6))
print("Trapezna formula daje:" , trapezna(lambda x : x*exp(-x)*cos(2*x),0,2*pi,6))
print("Simpsonova formula daje:" , Simpson(lambda x : x*exp(-x)*cos(2*x),0,2*pi,6))

Pravokutna formula daje: -0.06524222215454034
Trapezna formula daje: -0.2269939761611912
Simpsonova formula daje: -0.11915947349009065

Unutar Pythonovih biblioteka ima već ugrađenih funkcija za numeričku integraciju. Za detalje možete pogledati ovaj link

In [14]:
x = np.linspace(0,2*pi,7)
y=x*exp(-x)*cos(2*x)
print("Value  from Numpy.trapz:", np.trapz(y, x))

Value  from Numpy.trapz: -0.22699397616119124

Rekurzivno trapezno pravilo¶

Sada ćemo implementirati rekurzivno ptrapezno pravilo kako je to opisano u nastavnim materijalima.

In [15]:
def trap_rec(f,a,b, I0, k):
    if k==1:
        I1=(f(a)+f(b))*(b-a)/2.0
    else:
        n=2**(k-2)
        h=(b-a)/n
        x=a+h/2.0
        new_points_sum=0.0
        for j in range(n):
            new_points_sum=new_points_sum+f(x)
            x=x+h
        I1=(I0+h*new_points_sum)/2.0
    return I1
In [16]:
approx1=trap_rec(lambda x : x*exp(-x)*cos(2*x),0, 2*pi, 0.0,1)

print("Aproksimacija za k=1:", approx1)

approx2=trap_rec(lambda x : x*exp(-x)*cos(2*x),0, 2*pi,approx1,2)
print("Aproksimacija za k=2:", approx2)

approx3=trap_rec(lambda x : x*exp(-x)*cos(2*x),0, 2*pi,approx2,3)
print("Aproksimacija za k=3:", approx3)    

approx4=trap_rec(lambda x :x*exp(-x)*cos(2*x),0, 2*pi,approx3,4)

print("Aproksimacija za k=4:", approx4)  
approx5 =trap_rec(lambda x :x*exp(-x)*cos(2*x),0, 2*pi,approx4,5)
print("Aproksimacija za k=5:", approx5)    
    
approx6=trap_rec(lambda x :x*exp(-x)*cos(2*x),0, 2*pi,approx5,6)
print("Aproksimacija za k=6:", approx6)    
    
    

Aproksimacija za k=1: 0.03686184200729501
Aproksimacija za k=2: 0.44493519888808997
Aproksimacija za k=3: -0.35695084241677694
Aproksimacija za k=4: -0.17847542120838847
Aproksimacija za k=5: -0.13539684852415113
Aproksimacija za k=6: -0.12538580799307902

Implementirajmo funkciju koja računa rekurzivnim trapeznim pravilom vrijednost integrala do na zadanu točnost (pod time se podrazumijeva da je $|Ik-I_{k-1}|<\varepsilon$ gdje je $\varepsilon$ unaprijed dana točnost.

In [17]:
def trap_rec_eps(f,a,b, eps):
    I0=0.0
    I1= trap_rec(f,a,b,I0, 1)
    k=1
    while np.abs(I1-I0)>eps:
        I0=I1
        k=k+1
        I1= trap_rec(f,a,b,I0,k)
    return k, I1
In [10]:
k, approx=trap_rec_eps(lambda x :x*exp(-x)*cos(2*x),0, 2*pi,1.0e-8)
print("Aproksimacija je:", approx)
print("Broj nivoa (k) je:", k)

Aproksimacija je: -0.12212260771307237
Broj nivoa (k) je: 16

Rombergov algoritam¶

Sada ćemo primijeniti Richardsonovu ekstrapolaciju kako bismo poboljšali točnost kompozitnog trapeznog pravila. Takav postupak nazivamo Rombergov algoritam. Za primjenu Rombergovog algoritma možemo koristiti već ugrađenu funkciju

In [18]:
from scipy import integrate
from numpy import sin

result = integrate.romberg(sin, 0, pi, show=True)

Romberg integration of <function vectorize1.<locals>.vfunc at 0x000002824F924E18> from [0, 3.141592653589793]

 Steps  StepSize   Results
     1  3.141593  0.000000 
     2  1.570796  1.570796  2.094395 
     4  0.785398  1.896119  2.004560  1.998571 
     8  0.392699  1.974232  2.000269  1.999983  2.000006 
    16  0.196350  1.993570  2.000017  2.000000  2.000000  2.000000 
    32  0.098175  1.998393  2.000001  2.000000  2.000000  2.000000  2.000000 

The final result is 2.000000000001321 after 33 function evaluations.
In [ ]: