Numerical Integration Method
Soal
Soal Nomor 1
Diketahui: * pdf dari distribusi Kumaraswamy:
\[f(x|a,b) = ab x^{a-1}(1-x^a)^{b-1}\]
- a = 5
- b = 3
- batas 0 < x < 1
Substitusi nilai a dan b ke pdf \[f(x) = 5.3 x^{5-1}(1-x^5)^{3-1}\\ = 15 x^{4}(1-x^5)^{2}\]
Sehingga diperoleh E(x) \[E(x)=\int_0^1 x.f(x)\quad dx \\ =\int_0^1 x.15x^4(1-x^5)^2\quad dx\\=\int_0^1 15x^5(1-x^5)^2\quad dx\]
1.a Menghitung E(X) Secara Manual
Metode Trapezoida untuk n = 4
-Diketahui \[f(x) = 15 x^{5}(1-x^5)^{2},\quad a=0,\quad b=1\]
-Panjang sub interval \[h=\frac{b-a}{n}=\frac{1-0}{4}=\frac{1}{4}\] - Bentuk metode Trapezoida untuk n = 4
\[T_4 = \frac{h}{2}[f(x_0)+2f(x_1)+2f(x_2)+2f(x_3)+2f(x_4)]\] - Nilai f(x) untuk ‘xi’ \[ f(x_0)= f(0) = 0 \\f(x_1)=f(\frac{1}{4})= 0.01461984 \\f(x_2)=f(\frac{2}{4})= 0.4399109 \\f(x_3)=f(\frac{3}{4})= 2.070617 \\f(x_4)=f(\frac{4}{4})= 0\]
-Menghitung Nilai T_4 Karena h = 1, maka \[T_4 =\frac{1}{4}\frac{1}{2}[0+2(0.01461984)+2(0.4399109)+2(2.070617)+0)]=0.6312869\]
Sehingga \[ \int_0^1 15x^5(1-x^5)^2dx \quad \thickapprox \quad T_4 = 0.6312869\] Jadi, diperoleh E(x) menggunakan ‘metode Trapezoida’ untuk n = 4 yaitu E(x) ≈ 0.6312869
Metode Simpson untuk n = 4
-Diketahui \[f(x) = 15 x^{5}(1-x^5)^{2},\quad a=0,\quad b=1\]
-Panjang sub interval \[h=\frac{b-a}{n}=\frac{1-0}{4}=\frac{1}{4}\]
- Bentuk metode Simpson untuk n = 4
\[S_4 = \frac{h}{3}[f(x_0)+4f(x_1)+2f(x_2)+4f(x_3)+f(x_4)]\] - Nilai f(x) untuk ‘xi’ \[ f(x_0)= f(0) = 0 \\f(x_1)=f(\frac{1}{4})= 0.01461984 \\f(x_2)=f(\frac{2}{4})= 0.4399109 \\f(x_3)=f(\frac{3}{4})= 2.070617 \\f(x_4)=f(\frac{4}{4})= 0\]
-Menghitung Nilai S_4 Karena h = 1, maka \[S_4 =\frac{1}{4}\frac{1}{3}[0+4(0.01461984)+2(0.4399109)+4(2.070617)+0)]=0.7683974\]
Sehingga \[ \int_0^1 15x^5(1-x^5)^2dx \quad \thickapprox \quad S_4 = 0.7683974\] Jadi, diperoleh E(x) menggunakan ‘metode Simpson’ untuk n = 4 yaitu E(x) ≈ 0.7683974
Metode Four-Point Gauss Quadrature
Metode four-point gaussian quadratur, berarti nilai n=4. Berdasarkan tabel gaussian quadratur nilai-nilai koefisien dan titik gauss adalah sebagai berikut:
| n | Koefisien | Titik Gauss |
|---|---|---|
| 4 | C_1=0.3478548 | x_1=−0.861136311 |
| C_2=0.6521452 | x_2=−0.339981043 | |
| C_3=0.6521452 | x_3=0.339981043 | |
| 4 | C_4=0.3478548 | x_4=0.861136311 |
langkah 1 Fungsi f(x) dimana dari soal diminta menggunakan E(X) dan selang integrasinya (a,b)
telah diketahui yaitu:
\[f(x) = 15 x^5(1-x^5)^2,\quad a=0,\quad b=1\]
langkah 2 Sebelum kita menggunakan nilai-nilai koefisien dan titik gauss, perlu dilakukan transformasi sehingga domain integral menjadi [−1,1].
\[x=\frac{1}{2}[t(b-a)+(a+b)] \\=\frac{1}{2}[t(1-0)+(0+1)] \\=\frac{1}{2}(t+1) \\x=\frac{t+1}{2}\]
selanjutnya \[x=\frac{dx}{dt}=\frac{1}{2}\\ dx=\frac{1}{2}dt\] Substitusi nilai x ke persamaan \[ \int_0^1 15x^5(1-x^5)^2dx \quad =\int_{-1}^1 15\left[\frac{t+1}{2}\right]^5\left[1-\left[\frac{t+1}{2}\right]^5\right]^2.\frac{1}{2}dt \\=\int_{-1}^1 \frac{15}{2}\left[\frac{t+1}{2}\right]^5\left[1-\left[\frac{t+1}{2}\right]^5\right]^2dt\]
Menentukan nilai I \[\int_{-1}^1 \frac{15}{2}\left[\frac{t+1}{2}\right]^5\left[1-\left[\frac{t+1}{2}\right]^5\right]^2dt\quad \thickapprox \quad I \] \[I= C_1f(t_1)+C_2f(t_2)+C_3f(t_3)+C_4f(t_4) \\= 0.3478548.f(-0.86113631)+0.6521452.f(−0.339981043)\\+ 0.6521452.f(0.339981043)+0.3478548.f(0.86113631)\]
Menentukan Nilai f(x) untuk setiap x_i *\[f(-0.86113631)=\int_{-1}^1 \frac{15}{2}\left[\frac{-0.86113631+1}{2}\right]^5\left[1-\left[\frac{-0.86113631+1}{2}\right]^5\right]^2dt = 1.21x10^{-5}\]
*\[f(−0.339981043)=\int_{-1}^1\frac{15}{2}\left[\frac{−0.339981043+1}{2}\right]^5\left[1-\left[\frac{-0.339981043+1}{2}\right]^5\right]^2dt = 0.02912641\]
*\[f(0.339981043)=\int_{-1}^1\frac{15}{2}\left[\frac{0.339981043+1}{2}\right]^5\left[1-\left[\frac{0.339981043+1}{2}\right]^5\right]^2dt = 0.7575892\]
*\[f(0.86113631)=\int_{-1}^1\frac{15}{2}\left[\frac{0.86113631+1}{2}\right]^5\left[1-\left[\frac{0.86113631+1}{2}\right]^5\right]^2dt = 0.4779089\]
Sehingga \[I== 0.3478548(1.21x10^{-5})+0.6521452(0.02912641)\\+ 0.6521452(0.7575892)+0.3478548(0.4779089) \\=0.6792999\]
Jadi, \[\int_0^1 15x^5(1-x^5)^2dx\quad \thickapprox \quad I = 0.6792999\]
Komparasi Metode
Nilai_EX <- rbind(0.631286879,
0.768397357,
0.679299934
)Nama_Metode1 <- c("Metode Trapezoida untuk n=4",
"Metode Simpson untuk n=4",
"Metode four-point Gauss Quadrature"
)hasil <- data.frame(Nama_Metode1, Nilai_EX)
hasil## Nama_Metode1 Nilai_EX
## 1 Metode Trapezoida untuk n=4 0.6312869
## 2 Metode Simpson untuk n=4 0.7683974
## 3 Metode four-point Gauss Quadrature 0.67929991.b Menghitung E(X) dengan Nilai Exact
\[E(X)= \frac{bΓ(1+(\frac{1}{2})Γ(b))}{Γ(1+\frac{1}{a}+b)}\]
a <- 5
b <- 3
Eksak <- (b*gamma(1+(1/a))*gamma(b))/(gamma(1+(1/a)+b))
Eksak## [1] 0.7102273Fungsi
EX <- function(x){
x*a*b*x^(a-1)*(1-x^(a))^(b-1)
}# Menggambar kurva f(x)
curve(EX,0,1)
abline(v=seq(0,1,length.out = 7),col=77,lty=2)
Metode Trapezoida
mendefinisikan fungsi f(x)
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)
}f = function(x){15*x^5 - 30*x^10 + 15*x^15}
exact_value=0.7102273
tol <- 0.0001
err <- 1
n = 4
while(err>tol){
res_trap <- trapezoid(f,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-05Metode 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)
}f = function(x){15*x^5 - 30*x^10 + 15*x^15}
exact_value=0.7102273
tol <- 0.0001
err <- 1
n = 4
while(err>tol){
res_simp <- simpson_n(f,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-05Metode Gauss Quadrature
f = function(x){15*x^5 - 30*x^10 + 15*x^15}
exact_value=0.7102273
tol <- 0.0001
err <- 1
n = 4
library(pracma)
while(err>tol){
gL <- gaussLegendre(n = n,a = 0,1)
Ci <- gL$w # koefisien
xi <- gL$x # gauss point
res_gl <- sum(Ci * f(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-05Olehnya itu, untuk Metode Gauss Quadrature digunakan n=6
Komparasi Metode untuk tiap sub Interval Berikut perbandingan banyaknya sub interval n untuk mendapatkan E(x) dengan toleransi 0.0001.
exact <- 0.7102273
ertrap <- 0.7102273-0.6312869
ersimp <- 0.7683974-0.7102273
ergaquad <- 0.7102273-0.7101422
ErrorTollerance <- rbind(exact, ertrap,ersimp,ergaquad)
ErrorTollerance## [,1]
## exact 0.7102273
## ertrap 0.0789404
## ersimp 0.0581701
## ergaquad 0.0000851Nama_Metode <- c("Nilai Exact",
"Metode Trapezoida untuk n=4",
"Metode Simpson untuk n=4",
"Metode four-point Gauss Quadrature"
)library(tibble)
tibble(all_result <- data.frame(Nama_Metode, ErrorTollerance))## # A tibble: 4 x 2
## Nama_Metode ErrorTollerance
## <chr> <dbl>
## 1 Nilai Exact 0.710
## 2 Metode Trapezoida untuk n=4 0.0789
## 3 Metode Simpson untuk n=4 0.0582
## 4 Metode four-point Gauss Quadrature 0.0000851Dapat disimpulkan bahwa Metode terbaik yang digunakan adalah Metode four-points Gauss Quadrature dengan nilai error terkecil.
Soal Nomor 2
CDF dari suatu distribusi didefinisikan sebagai berikut:
\[F(x)=\int_0^\infty f(t)dt\]
Pdf dari distirbusi Eksponensial adalah \[f(x|λ)=λ e^{−λ}x\]
Menentukan CDF-nya dengan x = 4 dan λ=2 \[F(4)=\int_0^4f(x) dx \\f(t)=\int_0^42e^{-2}{t} dt\]
Menentukan Nilai Exact
a <- 0
b <- 4
F <- function(t){
2*exp(-2*t)
}
eksakcdf <- integrate(F, 0, 4)
eksakcdf## 0.9996645 with absolute error < 1.1e-14curve(F,0,4)
abline(v=seq(0,4,length.out = 7),col=7,lty=2)
Selanjutnya akan dicari cdf menggunakan 4 metode yaitu metode Trapezoida n=1000, Metode Simpson n = 300, Gauss Quadrature dengan n = 7, dan Monte Carlo dengan m = 10000. Untuk menghasilkan cdf akan digunakan nilai teloransi 0.00001
**1. Metode Trapezoida
trapezoid <- function(ftn, a, b, n = 1000) {
h <- (b-a)/n
h
x.vec <- seq(a, b, by = h)
f.vec <- sapply(x.vec, F) # ftn(x.vec)
Trap <- h*(f.vec[1]/2 + sum(f.vec[2:n]) + f.vec[n+1]/2)
return(Trap)
}
trapezoid(F,0,4,n = 1000)## [1] 0.9996699exact_value=0.9996645
tol <- 0.00001
err <- 1
n = 1000
while(err>tol){
res_trap <- trapezoid(F,0,4,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=1000, result=0.9996699, error=5.368911e-062. 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)
}
simpson_n(F,0,4,n = 300)## [1] 0.9996645exact_value=0.9996645
tol <- 0.0001
err <- 1
n = 500
while(err>tol){
res_simp <- simpson_n(F,0,4,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=500, result=0.9996645, error=3.773605e-083.Metode Gauss Quadrature
Menggunakan fungsi ‘gaussLegendre’ dari package ‘pracma’ R untuk n=7
gL <- gaussLegendre(n = 7,a = 0,4)
Ci <- gL$w # koefisien
xi <- gL$x # gauss point
I <- sum(Ci * F(xi))
I## [1] 0.9996645exact_value=0.9996645
tol <- 0.00001
err <- 1
n = 7
while(err>tol){
gL <- gaussLegendre(n = n,a = 0,4)
Ci <- gL$w # koefisien
xi <- gL$x # gauss point
res_gl <- sum(Ci * F(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=7, result=0.9996645, error=1.619311e-084. Metode Monte Carlo
x <- 4
lamda <- 2
F <- function(t){
lamda * exp(-lamda*t)
}# 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)
}set.seed(1001)
int_g <- mc_integral(F,0, 4,m = 5000)
int_g## [1] 1.001289Membandingkan Nilai Error
exact_value = 0.9996645
F = function(t){ 2*exp(-2*t) }
# Metode Trapezoida
result1 = trapezoid(F,a=0,b=4,n=1000)
err1 = abs(exact_value-result1)
# Metode Simpson
result2 = simpson_n(F,a=0,b=4,n=300)
err2 = abs(exact_value-result2)
# Metode Gauss Quadrature
result3 = sum(Ci * F(xi))
err3 = abs(exact_value-result3)
# Metode Monte Carlo
set.seed(1001)
result4=mc_integral(F,a=0,b=4,m=5000)
err4 = abs(exact_value-result4)
Metode = c("Trapezoida",
"Simpson",
"Gauss Quadrature",
"Monte Carlo")
n = c(1000,300,7,5000)
Hasil = c(result1,result2,result3,result4)
Error = c(err1,err2,err3,err4)komp.metode = data.frame(Metode,n, Hasil, Error)
komp.metode## Metode n Hasil Error
## 1 Trapezoida 1000 0.9996699 5.368911e-06
## 2 Simpson 300 0.9996645 4.018025e-08
## 3 Gauss Quadrature 7 0.9996645 1.619311e-08
## 4 Monte Carlo 5000 1.0012889 1.624354e-03Kesimpulan
Berdasarkan perhitungan fdf dari distribusi Eksponensial menggunakan 4 Metode yakni Metode Trapezoida n=1000, Metode Simpson n = 300, Gauss Quadrature dengan n = 7, dan Monte Carlo , dapat disimpulkan bahwa metode yang terbaik adalah metode Gauss Quadrature n=7, toleransi=0.00001 karena memiliki nilai error paling kecil dan hasil yang mendekati exactnya.