Metode Integral Numerik

Soal 1A

Diketahui pdf (probability density function) dari distribusi Kumaraswamy adalah sebagai berikut: \[𝑓(𝑥|𝑎,𝑏)=𝑎𝑏x^{𝑎−1}(1−𝑥^{𝑎})^{𝑏−1}, 𝑥∈(0,1)\]

Jika diketahui \(a=5\) dan \(b=3\). Hitunglah \(E(X)\) dengan menggunakan metode Trapezoida untuk \(n=4\), Simpson untuk \(n=4\) dan four-point Gauss Quadrature tanpa menggunakan software apapun. (Boleh diketik ataupun ditulis tangan kemudian difoto)

Metode Trapezoida untuk \(n=4\)

knitr::include_graphics("/Users/mac/Documents/Trapezoida.png")

Metode Simpson untuk \(n=4\)

knitr::include_graphics("/Users/mac/Documents/Simpson.png")

Metode four-point Gauss Quadrature

knitr::include_graphics("/Users/mac/Documents/Gaus Quadrature.png")

Perbandingan Perhitungan Excel

\(n\) Metode \(E(X)\)
4 Trapezoida 0.631287
4 Simpson 0.768397
4 Gauss Quadrature 0.679299

Soal 1B

Jika diketahui nilai exact dari \(𝐸(𝑋)=\frac {𝑏Γ(1+\frac{1}{a})Γ(𝑏)}{Γ(1+\frac{1}{a}+𝑏)}\) metode mana yang paling mendekati untuk menghitung nilai \(𝐸(𝑋)\) dengan toleransi 0.0001? Gunakan R!

Nilai exact dari \(𝐸(𝑋)=\frac {𝑏Γ(1+\frac{1}{a})Γ(𝑏)}{Γ(1+\frac{1}{a}+𝑏)} =\frac {3Γ(1+\frac{1}{5})Γ(3)}{Γ(1+\frac{1}{5}+3)}= 0.7102273\)

a <- 5
b <- 3
E_k <- (b*gamma(1+(1/a))*gamma(b))/(gamma(1+(1/a)+b))
E_k
## [1] 0.7102273

dengan toleransi 0.0001 nilai tersebut dibandingkan dengan tiga metode sebelumnya menggunakan software R

Metode Trapezoida

trapezoid <- function(ftn, a, b, n = 100) {
     h <- (b-a)/n
     x.vec <- seq(a, b, by = h)
     f.vec <- sapply(x.vec, ftn)     # ftn(x.vec)
     Trap <- h*(f.vec[1]/2 + sum(f.vec[2:n]) + f.vec[n+1]/2)
     return(Trap)
}
EXK <- function(x){
  x*a*b*x^(a-1)*(1-x^(a))^(b-1)
}

Menghitung integral menggunakan fungsi trapezoid

trapezoid(EXK,0,1,n = 4)
## [1] 0.6312869

Menggunakan fungsi trapzfun dari package pracma

trapzfun(EXK,0,1,maxit =  4)
## $value
## [1] 0.7098069
## 
## $iter
## [1] 4
## 
## $rel.err
## [1] 0.005901199

menggunakan R hingga hasil integral mendekati hasil exact-nya dengan toleransi 0.0001

exact_value=0.7102273
tol <-  0.0001
err <- 1
n = 4

while(err>tol){
  res_trap <- trapezoid(EXK,0,1,n = n)
  
  err <- abs(res_trap-exact_value)
  
  cat("n=",n,", result=",res_trap,", error=",err,"\n",sep = "")
  
  n=n+1
  if(n==1000){break}
}
## n=4, result=0.6312869, error=0.07894042
## n=5, result=0.6738123, error=0.03641497
## n=6, result=0.6914928, error=0.01873454
## n=7, result=0.6997126, error=0.01051473
## n=8, result=0.7039057, error=0.006321603
## n=9, result=0.7062116, error=0.004015656
## n=10, result=0.7075597, error=0.002667593
## n=11, result=0.7083885, error=0.001838822
## n=12, result=0.7089199, error=0.001307439
## n=13, result=0.7092729, error=0.0009544148
## n=14, result=0.7095147, error=0.000712649
## n=15, result=0.7096846, error=0.0005426739
## n=16, result=0.7098069, error=0.0004204042
## n=17, result=0.7098966, error=0.0003306616
## n=18, result=0.7099637, error=0.0002636072
## n=19, result=0.7100146, error=0.0002127014
## n=20, result=0.7100538, error=0.0001734994
## n=21, result=0.7100844, error=0.0001429188
## n=22, result=0.7101085, error=0.0001187833
## n=23, result=0.7101278, error=9.953088e-05

Metode Simpson

simpson_n <- function(ftn, a, b, n = 100) {
    n <- max(c(2*(n %/% 2), 4))
    h <- (b-a)/n
    x.vec1 <- seq(a+h, b-h, by = 2*h)       # ganjil
    x.vec2 <- seq(a+2*h, b-2*h, by = 2*h)   # genap
    f.vec1 <- sapply(x.vec1, ftn)   # ganjil
    f.vec2 <- sapply(x.vec2, ftn)   # genap
    S <- h/3*(ftn(a) + ftn(b) + 4*sum(f.vec1) + 2*sum(f.vec2))
    return(S)
}

Menghitung integral menggunakan fungsi Simpson

simpson_n(EXK,0,1,n = 4)
## [1] 0.7683974

menggunakan R hingga hasil integral mendekati hasil exact-nya dengan toleransi 0.0001

exact_value=0.7102273
tol <-  0.0001
err <- 1
n = 4

while(err>tol){
  res_simp <- simpson_n(EXK,0,1,n = n)
  
  err <- abs(res_simp-exact_value)
  
  cat("n=",n,", result=",res_simp,", error=",err,"\n",sep = "")
  
  n=n+1
  if(n==1000){
    break
  }
  
}
## n=4, result=0.7683974, error=0.05817006
## n=5, result=0.7683974, error=0.05817006
## n=6, result=0.7497082, error=0.03948094
## n=7, result=0.7497082, error=0.03948094
## n=8, result=0.728112, error=0.01788467
## n=9, result=0.728112, error=0.01788467
## n=10, result=0.7188088, error=0.008581532
## n=11, result=0.7188088, error=0.008581532
## n=12, result=0.7147289, error=0.004501593
## n=13, result=0.7147289, error=0.004501593
## n=14, result=0.712782, error=0.002554712
## n=15, result=0.712782, error=0.002554712
## n=16, result=0.711774, error=0.001546662
## n=17, result=0.711774, error=0.001546662
## n=18, result=0.7112144, error=0.0009870756
## n=19, result=0.7112144, error=0.0009870756
## n=20, result=0.7108852, error=0.0006578651
## n=21, result=0.7108852, error=0.0006578651
## n=22, result=0.7106819, error=0.0004545628
## n=23, result=0.7106819, error=0.0004545628
## n=24, result=0.7105511, error=0.0003237819
## n=25, result=0.7105511, error=0.0003237819
## n=26, result=0.710464, error=0.0002366803
## n=27, result=0.710464, error=0.0002366803
## n=28, result=0.7104042, error=0.0001769138
## n=29, result=0.7104042, error=0.0001769138
## n=30, result=0.7103621, error=0.0001348299
## n=31, result=0.7103621, error=0.0001348299
## n=32, result=0.7103318, error=0.0001045198
## n=33, result=0.7103318, error=0.0001045198
## n=34, result=0.7103096, error=8.225045e-05

Metode four-point Gauss Quadrature

Menggunakan fungsi gaussLegendre dengan domain asal

Legendre order 4

gL <- gaussLegendre(n = 4,a = 0,1)
gL
## $x
## [1] 0.06943184 0.33000948 0.66999052 0.93056816
## 
## $w
## [1] 0.1739274 0.3260726 0.3260726 0.1739274

Mendefinisikan koefisien dan gauss point lalu menghitung integral

Ci <- gL$w # koefisien
Ci
## [1] 0.1739274 0.3260726 0.3260726 0.1739274
xi <- gL$x # gauss point
xi
## [1] 0.06943184 0.33000948 0.66999052 0.93056816
I <- sum(Ci * EXK(xi))
I
## [1] 0.6792999

menggunakan R hingga hasil integral mendekati hasil exact-nya dengan toleransi 0.0001

exact_value=0.7102273
tol <-  0.0001
err <- 1
n = 4

while(err>tol){
  
  gL <- gaussLegendre(n = n,a = 0,1)

  Ci <- gL$w # koefisien
  xi <- gL$x # gauss point
  
  res_gl <- sum(Ci * EXK(xi))

  err <- abs(res_gl-exact_value)
  
  cat("n=",n,", result=",res_gl,", error=",err,"\n",sep = "")
  
  n=n+1
  if(n==1000){
    break
  }
  
}
## n=4, result=0.6792999, error=0.03092738
## n=5, result=0.707375, error=0.002852279
## n=6, result=0.7101422, error=8.507726e-05

Perbandingan Perhitungan R

\(n\) Metode Hasil Error
23 Trapezoida 0.7101278 9.950361e-05
34 Simpson 0.7103096 8.227772e-05
6 Gauss Quadrature 0.7101422 8.507726e-05

Nilai exact dari \(E(X)\) sebesar 0.7102273.

Jika dilihat dari tabel di atas, maka hasil yang paling mendekati nilai exact \(𝐸(𝑋)\) dengan toleransi 0.0001 adalah metode Simpson karena errornya paling kecil.

Soal 2

Diketahui pdf (probability density function) dari distribusi Eksponensial adalah sebagai berikut: \[ 𝑓(𝑥|𝑎,𝑏)=𝜆 e^{−𝜆𝑥}, 𝑥∈(0,∞) \]

Jika diketahui \(𝜆=2\). Hitunglah CDF dari distribusi exponensial tersebut untuk \(x=4\), dengan menggunakan metode Trapezoida, Simpson, Gauss Quadrature dan Monte Carlo dengan R. Metode mana yang paling baik? Silahkan tentukan nilai toleransi, \(n\) dan \(m\) sendiri.

Fungsi distribusi kumulatif (CDF) distribusi Eksponensial untuk \(x=4\), jika diketahui \(λ=2\) adalah:

\[ F(x=4|λ=2)=\int_{0}^{4} 2e^{-2t} \] dengan nilai exact:

\[ F(x)=1-e^{-λt}= 1-e^{-(2)(4)} = 0.9996645\]

a <- 0
b <- 4
F <- function(t){
  2*exp(-2*t)
}
F_int <- integrate(F, 0, 4)
F_int
## 0.9996645 with absolute error < 1.1e-14

atau

cdf = 1-exp(-8)
cdf
## [1] 0.9996645

Fungsi untuk Metode Monte Carlo

mc_integral <- function(ftn, a, b, m=1000){
  #Membangkitkan x berdistribusi U(a,b)
  x <- runif(m,a,b)
  # Menghitung rata-rata dari output fungsi
  Gx <- ftn(x)
  Gx_m <- mean(Gx)
  theta.hat <- (b-a)*Gx_m
  return(theta.hat)
}

Menghitung CDF

Penghitungan CDF digunakan nilai \(n=250\) untuk metode Trapezoid dan Simpson, \(n=8\) untuk Gauss Quadrature, contoh acak \(m=1000\), dan toleransi 0.0001

cdf = 1-exp(-8)
F = function(t){2*exp(-2*t) }

# Metode Trapezoid
ht  = trapezoid(F,0,4,n=250)
er1 = abs(cdf-ht)

# Metode Simpson
hs  = simpson_n(F,0,4,n=250)
er2 = abs(cdf-hs)

# Metode Gauss Quadrature
gL  = gaussLegendre(n = 8,a = 0,4)
Ci  = gL$w # koefisien
xi  = gL$x # gauss point
hgq  = sum(Ci * F(xi))
er3 = abs(hgq-cdf)

# Metode Monte Carlo
set.seed(123)
m=1000
hmc=mc_integral(F,a=0,b=4,m=m)
er4 = abs(hmc-cdf)

Metode = c("Trapezoid","Simpson","Gauss Quadrature","Monte Carlo")
n      = c(250,250,8,m)
Hasil  = c(ht,hs,hgq,hmc)
Error  = c(er1,er2,er3,er4)

data.frame(Metode,n,Hasil,Error)
##             Metode    n     Hasil        Error
## 1        Trapezoid  250 0.9997498 8.530325e-05
## 2          Simpson  250 0.9996645 5.822758e-09
## 3 Gauss Quadrature    8 0.9996645 8.724643e-10
## 4      Monte Carlo 1000 1.0018262 2.161698e-03

Jika dilihat dari tabel di atas, maka hasil yang paling mendekati nilai exact \(F(𝑋) = 0.9996645\) dengan toleransi 0.0001 adalah metode Gaus Quadratur karena errornya paling kecil.