Numerička integracija (nastavak)¶

U prethodnom je predavanju postavljen video s rješenjima zadataka. Pogledajte ga!

U sljedećem dijelu ćemo pokušati kompozitnu trapeznu formulu zapisati rekurzivno, kako bismo imali pšto manje izvrjednjavanja funkcije kod povećanja broja točaka!

In [25]:
from IPython.lib.display import YouTubeVideo
vid = YouTubeVideo("n5JAROOhSJA")
display(vid)

Slijedi usporedba sa vrijednošću koju dobijemo korištenjem ugrađene funkcije iz scipy biblioteke

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

import numpy as np
from numpy import exp
from numpy import cos
from numpy import pi
from numpy import sin

#trapez(lambda x : x*exp(-x)*cos(2*x),0,2*pi,5)

approx1=trapez(lambda x : sin(x),0, pi,8)
print("Aproksimacija-implementacija=", approx1)

AProksimacija= 1.9742316019455508
In [23]:
from scipy import integrate
x = np.linspace(0,pi,9)
y=sin(x)
print("Value  from scipy(trapezoid):", integrate.trapz(y, x))

Value  from scipy(trapezoid): 1.9742316019455508

U sljedećem kodu implementira se rekurzivno trapezno pravilo i rješavamo primjer koji je riješen u video predavanjima

In [18]:
#funkcija koja računa vrojednost integrala I_k sa brojem intervala 2^{k-1} iz vrijednosti integrala I_{k-1}


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
# Primjer s video predavanja

approx2=trap_rec(lambda x : sin(x),0, pi,approx1,5)
print("Aproksimacija za k=5:", approx2)

approx3=trap_rec(lambda x : sin(x),0, pi,approx2,6)
print("Aproksimacija za k=6:", approx3)    

approx4=trap_rec(lambda x : sin(x),0, pi,approx3,7)
print("Aproksimacija za k=7:", approx4)    
    

Aproksimacija za k=5: 1.9935703437723393
Aproksimacija za k=6: 1.998393360970144
Aproksimacija za k=7: 1.9995983886400366

D.Z. implementirati funkciju koja računa rekurzivnim trapeznim pravilo 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. Ako se ne može postići navedeno, ispisati poruku da postupak ne može dati takvu točnost)

Richardsonova ekstrapolacija¶

Richardsonova ekstrapolacija općenit je postupak poboljšavanja točnosti aproksimacijskih formula. U sljedećem videu bit će objašnjen općenit postupak i primjena na numeričko deriviranje

In [26]:
vid = YouTubeVideo("y2GeEeY54EY")
display(vid)

Rombergov algoritam¶

Sada ćemo primijeniti richardsonovu ekstrapolaciju kako bismo poboljšali točnost kompozitnog trapeznog pravila. Takav postupak nazivamo Rombergov algoritam.

In [27]:
vid = YouTubeVideo("r0B5vZopCKw")
display(vid)

Za dz. razmisliti o implementaciji Rombergovog algoritma, ali može se koristiti i već ugrađena funkcija

In [21]:
from scipy import integrate

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

Romberg integration of <function vectorize1.<locals>.vfunc at 0x00000159FCFC56A8> 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 [ ]: