Tugas Praktikum 10
library(pracma)Soal 1
1a
E(x) dengan metode Trapezoida
E(x) dengan metode Simpson
E(x) dengan metode Gauss-Quadrature
Perbandingan Hasil perhitungan Manual tiga Metode dengan n=4.
1b
Pada soal 1b) akan dihitung nilai E(X) menggunakan metode trapezoida, Simpson dan Gauss Quadrature dengan nilai toleransi sebesar 0.0001
a<-5
b<-3
meang <- b * gamma (1+1/a) * gamma (b) / gamma (1+1/a+b)
meang## [1] 0.7102273Maka Exact Value = 0.7102273
Metode Trapezoid
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
ft <- function(x){
15*x^5 - 30*x^10 + 15*x^15
}
gL <- gaussLegendre(n = 4,a = 0,1)
Ci <- gL$w # koefisien
Ci## [1] 0.1739274 0.3260726 0.3260726 0.1739274xi <- gL$x # gauss point
xi## [1] 0.06943184 0.33000948 0.66999052 0.93056816I <- sum(Ci * ft(xi))
I## [1] 0.6792999exact_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 * ft(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-05Komparasi Hasil perhitungan tiga metode
Perbandingan banyaknya sub interval n untuk mendapatkan E(x) dengan nilai toleransi 0.0001
data.frame(Metode=c("Trapezoida","Simpson","Gauss Quadrature"),"Interval"=c(23,34,6))## Metode Interval
## 1 Trapezoida 23
## 2 Simpson 34
## 3 Gauss Quadrature 6Berdasarkan hasil iterasi metode yang terbaik adalah metode Gauss Quadrature karena hanya membutuhkan 6 sub interval.
Soal 2
a <- 0
b <- 4
lambda <- 2
fall <- function(x, lambda) lambda*exp(-lambda*x)
integrate(fall, a, b, lambda)## 0.9996645 with absolute error < 1.1e-14Nilai eksak = 0.9996645
Dengan nilai exact =0.9996645, Menggunakan metode Trapezoida dan Simpson untuk n = 100, Gauss Quadrature dengan n = 6, dan Monte Carlo dengan m = 1000 akan diperoleh hasil yang akan dikomparasi.
Fungsi Integral Monte Carlo dengan m=1000
mc_integral <- function(fall, a, b,m=1000){
#Membangkitkan x berdistribusi U(a,b)
x <- runif(m,a,b)
# Menghitung rata-rata dari output fungsi
Gx <- fall(x)
Gx_m <- mean(Gx)
theta.hat <- (b-a)*Gx_m
return(theta.hat)
}Menghitung CDF dari distribusi exponensial.
exact_value = 0.9996645
fp = function(x){ 2*exp(-2*x) }
# Metode Trapezoida
ht = trapezoid(fp,a=0,b=4,n=100)
e1 = abs(exact_value -ht)
# Metode Simpson
hs = simpson_n(fp,a=0,b=4,n=100)
e2 = abs(exact_value -hs)
# Metode Gauss Quadrature
gL = gaussLegendre(n = 100,a = 0,4)
Ci = gL$w # koefisien
xi = gL$x # gauss point
hg = sum(Ci * fp(xi))
e3 = abs(hg-exact_value )
# Metode Monte Carlo
set.seed(41)
m=1000
hmc=mc_integral(fp,a=0,b=4,m=m)
e4 = abs(hmc-exact_value )
Metode = c("Trapezoida","Simpson","Gauss Quadrature","Monte Carlo")
n = c(100,100,100,m)
Hasil = c(ht,hs,hg,hmc)
Error = c(e1,e2,e3,e4)
data.frame(Metode,n,Hasil,Error)## Metode n Hasil Error
## 1 Trapezoida 100 1.0001976 5.331349e-04
## 2 Simpson 100 0.9996648 2.646781e-07
## 3 Gauss Quadrature 100 0.9996645 3.737209e-08
## 4 Monte Carlo 1000 0.9539666 4.569792e-02Berdasarkan komparasi dari ke-empat metode tersebut, metode yang terbaik adalah gauss quadrature dengan hasil paling mendekati exact value dengan nilai error terkecil.