Numeryczne metody obliczania całek

Wstęp

Proste metody całkowania numerycznego polegają na przybliżeniu całki za pomocą odpowiedniej sumy ważonej wartości całkowanej funkcji w kilku punktach. Aby uzyskać dokładniejsze przybliżenie dzieli się przedział całkowania na niewielkie fragmenty. Ostateczny wynik jest sumą oszacowań całek w poszczególnych podprzedziałach. Najczęściej przedział dzieli się na równe podprzedziały, ale bardziej wyszukane algorytmy potrafią dostosowywać krok do szybkości zmienności funkcji.

W niniejszej pracy zostaną przedstawione trzy moetody numeryczne:

• metoda prostokątów
• metoda trapezów
• metoda Monte Carlo

Stworzone algorytmy do obliczania całek, przetestujemy na dwóch przykładowych funkcjacjach:

eq1=function(x){2*x^3-2*x^2+x-2}
eq2=function(x){sin(x)}

Dla obu funkcji, zakres całkowania ustalmy w przedziale x <-2,3>, zatem wartości całek, powinny odpowiednio wynosić:

integrate(eq1,-2,3)

## 1.666667 with absolute error < 4.9e-13

oraz

integrate(eq2,-2,3)

## 0.5738457 with absolute error < 3.8e-14

Metoda prostokątów

W metodzie prostokątów (ang. rectangular integration) korzystamy z definicji całki oznaczonej Riemanna , w której wartość całki interpretowana jest jako suma pól obszarów pod wykresem krzywej w zadanym przedziale całkowania . Sumę tę przybliżamy przy pomocy sumy pól odpowiednio dobranych prostokątów.

Implementacja metody w programie RStudio:

midpoint=function(fun,x1,x2,n){
    dx=(x2-x1)/(n)
    dx2=dx/2
    x=seq(x1+dx2,x2-dx2,length=n)
    
     
    curve(fun(x),xlim=c(x1-0.1*n*dx,x2+0.1*n*dx), col="red", lwd=3)
    y=0 
    
    for(i in 1:(n)){
         y=y+fun(x[i])
         cord=c(x[i]-dx2,x[i]+dx2)
         polygon(c(cord,rev(cord)),c(fun(x[i]),fun(x[i]),0,0),col="gray")
    }
     
    return (y*dx)
    }

Przykład 1

midpoint(eq1,-2,3,50)

## [1] 1.6625

Przykład 2

midpoint(eq2,-2,3,20)

## [1] 0.5753428

Metoda trapezów

Opisana w poprzednim rozdziale metoda prostokątów nie jest zbyt dokładna, ponieważ pola użytych w niej prostokątów źle odwzorowują pole pod krzywą. Dużo lepszym rozwiązaniem jest zastosowanie zamiast nich trapezów o wysokości dx i podstawach równych odpowiednio wartości funkcji w punktach krańcowych.. Sama zasada nie zmienia się. Przybliżona wartość całki jest sumą pól wszystkich otrzymanych w ten sposób trapezów.

Implementacja metody w programie RStudio:

trapez=function(fun,x1,x2,n){
    x=seq(x1,x2,length=n)
    dx=(x2-x1)/(n-1)
     
    curve(fun(x),xlim=c(x1-0.1*n*dx,x2+0.1*n*dx), col="red", lwd=3)
     
    y=-fun(x1)/2
     
    for(i in 1:(n-1)){
         y=y+fun(x[i])
         cord=c(x[i],x[i+1])
         polygon(c(cord,rev(cord)),c(fun(cord),0,0),col="gray")
    }
     
    y=y+fun(x2)/2
           
    return (y*dx)
    }

Przykład 1

trapez(eq1,-2,3,50)

## [1] 1.675344

Przykład 2

trapez(eq2,-2,3,20)

## [1] 0.5705302

Metoda Monte Carlo

Do przybliżonego obliczania całki oznaczonej można również wykorzystać metody probabilistyczne. Należy pamiętać jednak, że wynik takiego całkowania jest też zmienną losową.

Metoda Monte Carlo (MC) jest stosowana do modelowania matematycznego procesów zbyt złożonych (obliczania całek, łańcuchów procesów statystycznych), aby można było przewidzieć ich wyniki za pomocą podejścia analitycznego. Istotną rolę w metodzie MC odgrywa losowanie (wybór przypadkowy) wielkości charakteryzujących proces, przy czym losowanie dokonywane jest zgodnie z rozkładem, który musi być znany.

Metoda została opracowana i pierwszy raz zastosowana przez Stanisława Ulama.

Implementacja metody w programie RStudio:

monte=function(fun,x1,x2,n){
    x=seq(x1,x2,length=n)
    y=fun(x)
    
    ymin=min(y)
    ymax=max(y)
    
    plot(x,y,type="l",col="blue", xlim=c(x1,x2), ylim=c(ymin,ymax),lwd=3)
    
    py=runif(n,ymin,ymax)
    px=runif(n,x1,x2)
    
    y=fun(px)
    z1=(y>0 & py>0 & py<=y)
    z2=(y<0 & py<0 & py>=y)
    
    points(px[z1],py[z1],col="green",pch='.')   
    points(px[z2],py[z2],col="red",pch='.')    
    points(px[!z1&!z2],py[!z1&!z2],col="black",pch='.')
    
    
    return(abs(x1-x2)*abs(ymax-ymin)*(sum(z1)-sum(z2))/n)
}

Przykład 1

monte(eq1,-2,3,75000)

## [1] 1.178667

Przykład 2

monte(eq2,-2,3,75000)

## [1] 0.5444