Integrasi Numerik

Titin Agustina ¹

2021-11-16

Soal

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)

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!

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.

Jawaban 1a

Metode Trapezoida untuk \(n=4\)

Masukkan nilai \(a\) dan \(b\) ke dalam pdf distribusi Kumaraswamy, sehingga diperoleh:

\[ 𝑓(𝑥|5,3)=15x^{4}(1−𝑥^{5})^{2} \]

Nilai \(E(X)\) dari \([0,1]\) adalah: \[E(X) = \int_{0}^{1} x \cdot f(x)dx = \int_{0}^{1} 15x^{5}(1−𝑥^{5})^{2} dx \]

Panjang setiap sub interval adalah sebagai berikut:

\[h=\frac{b-a}{n}=\frac{1-0}{4}=\frac{1}{4}\] Karena \(n=4\) maka bentuk dari metode trapezoidal adalah:

\[T_{4} = \frac{h}{2} [E(X_{0})+2E(X_{1})+2E(X_{2})+2E(X_{3})+E(X_{4})]\]

Menghitung fungsi \(E(X)\) pada setiap \(X_{i}\):

\(E(X_{0}) = E(0) = 15 (0)^{5}(1−(0)^{5})^{2} =0\)

\(E(X_{1}) = E(\frac{1}{4}) = 15 (\frac{1}{4})^{5}(1−(\frac{1}{4})^{5})^{2} =0.014619841\)

\(E(X_{2}) = E(\frac{1}{2}) = 15 (\frac{1}{2})^{5}(1−(\frac{1}{2})^{5})^{2} =0.439910889\)

\(E(X_{3}) = E(\frac{3}{4}) = 15 (\frac{3}{4})^{5}(1−(\frac{3}{4})^{5})^{2} =2.0706167866\)

\(E(X_{4}) = E(1) = 15 (1)^{5}(1−(1)^{5})^{2} =0\)

karena \(h=\frac{1}{4}\), maka diperoleh:

\[T_{4} = \frac{(\frac{1}{4})}{2} [0+2 (0.014619841)+2 (0.439910889)+2 (2.0706167866)+0] = 0.631286879\]

Sehingga

\[E(X) = \int_{0}^{1} 15x^{5}(1−𝑥^{5})^{2} dx \approx T_{4} = 0.631286879\]

Metode Simpson untuk \(n=4\)

Bentuk dari metode Simpson adalah:

\[S_{4} = \frac{h}{3} [E(X_{0})+4E(X_{1})+2E(X_{2})+4E(X_{3})+E(X_{4})]\]

karena \(h=\frac{1}{4}\), maka diperoleh

\[S_{4} = \frac{(\frac{1}{4})}{3} [0 + 4 (0.014619841) + 2 (0.439910889)+ 4 (2.070616786) + 0] = 0.768397357\]

Sehingga

\[ E(X) = \int_{0}^{1} 15x^{5}(1−𝑥^{5})^{2} dx \approx S_{4} = 0.7684 \]

Metode four-point Gauss Quadrature

Bentuk umum dari Gauss Quadrature adalah:

\[\int_{a}^{b}f(x)dx \approx I= \sum_{i=1}^{n}C_{i}f(x_{i})\]

dimana koefisien \(C_{i}\) merupakan bobot dan \(x_{i}\) merupakan titik Gauss yang berada dalam interval \([-1,1]\).

Nilai bobot \(C_{i}\) dan titik Gauss \(x_{i}\) dapat diperoleh saat domain integralnya adalah \([-1,1]\) dan \(f(x)\) merupakan polynomial berderajat \(n\). Untuk \(n=4\), maka nilai-nilai tersebut bisa dilihat pada tabel berikut:

\(n\)	Koefisien \(C_i\)	Titik Gauss \(x_i\)
4	\(C_1\) = 0.3478548	\(x_1\) = -0.86113631
	\(C_2\) = 0.6521452	\(x_2\) = -0.33998104
	\(C_3\) = 0.6521452	\(x_3\) = 0.33998104
	\(C_4\) = 0.3478548	\(x_4\) = 0.86113631

Jika domain integral berupa \([a,b]\), maka perlu ada transformasi sedemikian rupa sehingga domain integralnya menjadi \([-1,1]\)

\[\int_{a}^{b}f(x)dx \space \rightarrow \space \int_{-1}^{1}f(t)dt \]

Transformasi ini dilakukan dengan mengubah variabel sebagai berikut:

\[x=\frac{1}{2}[t(b-a)+a+b]\]

dan

\[dx=\frac{1}{2}(b-a)dt\]

Dengan menggunakan \(f(x)\), \(a\) dan \(b\) yang sama seperti pada dua metode sebelumnya, maka diperoleh:

\[x=\frac{1}{2}[t(1-0)+0+1]=\frac{t+1}{2}\] dan \[dx=\frac{1}{2}(1-0)dt=\frac{1}{2}dt\]

Sehingga integralnya menjadi:

\[ E(X) = \int_{0}^{1} 15x^{5}(1−𝑥^{5})^{2} dx = \int_{-1}^{1} 15(\frac{t+1}{2})^{5}(1−(\frac{t+1}{2})^{5})^{2} \frac{1}{2}dt \approx I \]

dimana

\(I = C_1E(t_1)+C_2E(t_2)+C_3E(t_3)+C_4E(t_4)\)

\(I = 0.3478548 \space E(−0.861136311) + 0.6521452 \space E(−0.339981043) + 0.6521452 \space E(0.339981043) + 0.3478548 \space E(0.861136311)\)

Setelah itu, masukkan nilai-nilai koefisien dan titik gauss saat \(n=4\) ke fungsi \(E(X)\) pada setiap \(x_i\)

\(E(-0.861136311) = \frac{15}{2} \left(\frac{-0.861136311 + 1}{2}\right)^{5}\left(1−\left(\frac{-0.861136311 + 1}{2}\right)^{5}\right)^{2} = 0.0000121019\)

\(E(−0.339981043) = \frac{15}{2} \left(\frac{−0.339981043 + 1}{2}\right)^{5}\left(1−\left(\frac{−0.339981043 + 1}{2}\right)^{5}\right)^{2} = 0.029126407\)

\(E(0.339981043) = \frac{15}{2} \left(\frac{0.339981043 + 1}{2}\right)^{5}\left(1−\left(\frac{0.339981043 + 1}{2}\right)^{5}\right)^{2} = 0.757589241\)

\(E(0.861136311) = \frac{15}{2} \left(\frac{0.861136311 + 1}{2}\right)^{5}\left(1−\left(\frac{0.861136311 + 1}{2}\right)^{5}\right)^{2} = 0.477908858\)

sehingga

\(I = 0.3478548 \space (0.0000121019) + 0.6521452 \space (0.029126407) + 0.6521452 \space (0.757589241) + 0.3478548 \space (0.477908858)\)

\(I = 0.679299934\)

Dengan demikian diperoleh:

\[ E(X) = \int_{0}^{1} 15x^{5}(1−𝑥^{5})^{2} dx = \int_{-1}^{1} 15(\frac{t+1}{2})^{5}(1−(\frac{t+1}{2})^{5})^{2} \frac{1}{2}dt \approx I = 0.679299934 \]

Tabel Hasil Perhitungan Manual

\(n\)	Metode	\(E(X)\)
4	Trapezoida	0.631286879
4	Simpson	0.768397357
4	Gauss Quadrature	0.679299934

Jawaban 1b

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

Tabel Hasil 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.

Jawaban 2

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=100\) untuk metode Trapezoid dan Simpson, \(n=7\) untuk Gauss Quadrature, contoh acak \(m=10000\), dan toleransi 0.0001

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

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

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

# Metode Gauss Quadrature
gL  = gaussLegendre(n = 7,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(308)
m=10000
hmc=mc_integral(F,a=0,b=4,m=m)
er4 = abs(hmc-cdf)

Metode = c("Trapezoid","Simpson","Gauss Quadrature","Monte Carlo")
n      = c(100,100,7,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   100 1.0001976 5.330976e-04
## 2          Simpson   100 0.9996648 2.273060e-07
## 3 Gauss Quadrature     7 0.9996645 5.356521e-08
## 4      Monte Carlo 10000 1.0284483 2.878372e-02

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

Terima kasih

Titin Agustina | G1501211005↩︎