P10 STA561 Metode Integral Numerik

SOAL 1

Alt text

Bagian a

Metode Trapezoida

Alt text

Metode Simpson

Alt text

Metode Gauss Quadrature

Alt text

Berikut adalah tabel gaussian quadrature nilai-nilai koefisien dan titik gauss dengan n = 4

Alt text

Bagian b

Dengan nilai toleransi = 0.0001, maka didapatkan nilai E(x) secara Exact sebagai berikut :

a<-5
b<-3
value <- b * gamma (1+1/a) * gamma (b) / gamma (1+1/a+b)
value

## [1] 0.7102273

Metode Trapezoida

# fungsi 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)
}

# mendefinisikan fungsi
E_mt <- function(x){
  (15*(x^5))*((1-(x^5))^2)
}

# menghitung nilai ekspetasi
exact_value=0.710227
tol <-  0.0001
err <- 1
n = 4

while(err>tol){
  res_trap <- trapezoid(E_mt,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.07894012
## n=5, result=0.6738123, error=0.03641467
## n=6, result=0.6914928, error=0.01873424
## n=7, result=0.6997126, error=0.01051443
## n=8, result=0.7039057, error=0.006321303
## n=9, result=0.7062116, error=0.004015356
## n=10, result=0.7075597, error=0.002667293
## n=11, result=0.7083885, error=0.001838522
## n=12, result=0.7089199, error=0.001307139
## n=13, result=0.7092729, error=0.0009541148
## n=14, result=0.7095147, error=0.000712349
## n=15, result=0.7096846, error=0.0005423739
## n=16, result=0.7098069, error=0.0004201042
## n=17, result=0.7098966, error=0.0003303616
## n=18, result=0.7099637, error=0.0002633072
## n=19, result=0.7100146, error=0.0002124014
## n=20, result=0.7100538, error=0.0001731994
## n=21, result=0.7100844, error=0.0001426188
## n=22, result=0.7101085, error=0.0001184833
## n=23, result=0.7101278, error=9.923088e-05

Didapatkan nilai E(x) sebesar 0.7101278 pada iterasi ke-23

Metode Simpson

# fungsi 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)
}

# mendefinisikan fungsi
E_ms <- function(x){
  (15*(x^5))*((1-(x^5))^2)
}


# menghitung nilai ekspetasi
exact_value=0.710227
tol <-  0.0001
err <- 1
n = 4

while(err>tol){
  res_simp <- simpson_n(E_ms ,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.05817036
## n=5, result=0.7683974, error=0.05817036
## n=6, result=0.7497082, error=0.03948124
## n=7, result=0.7497082, error=0.03948124
## n=8, result=0.728112, error=0.01788497
## n=9, result=0.728112, error=0.01788497
## n=10, result=0.7188088, error=0.008581832
## n=11, result=0.7188088, error=0.008581832
## n=12, result=0.7147289, error=0.004501893
## n=13, result=0.7147289, error=0.004501893
## n=14, result=0.712782, error=0.002555012
## n=15, result=0.712782, error=0.002555012
## n=16, result=0.711774, error=0.001546962
## n=17, result=0.711774, error=0.001546962
## n=18, result=0.7112144, error=0.0009873756
## n=19, result=0.7112144, error=0.0009873756
## n=20, result=0.7108852, error=0.0006581651
## n=21, result=0.7108852, error=0.0006581651
## n=22, result=0.7106819, error=0.0004548628
## n=23, result=0.7106819, error=0.0004548628
## n=24, result=0.7105511, error=0.0003240819
## n=25, result=0.7105511, error=0.0003240819
## n=26, result=0.710464, error=0.0002369803
## n=27, result=0.710464, error=0.0002369803
## n=28, result=0.7104042, error=0.0001772138
## n=29, result=0.7104042, error=0.0001772138
## n=30, result=0.7103621, error=0.0001351299
## n=31, result=0.7103621, error=0.0001351299
## n=32, result=0.7103318, error=0.0001048198
## n=33, result=0.7103318, error=0.0001048198
## n=34, result=0.7103096, error=8.255045e-05

Didapatkan nilai E(x) sebesar 0.7103096 pada iterasi ke-34

Metode Gauss Quadrature

library(pracma)

E_mfpgq <- function(x){
  (15*(x^5))*((1-(x^5))^2)
}

exact_value=0.710227
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 * E_mfpgq(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.03092708
## n=5, result=0.707375, error=0.002851979
## n=6, result=0.7101422, error=8.477726e-05

Didapatkan nilai E(x) sebesar 0.7101422 pada iterasi ke-6

Perbandingan Metode

Alt text

Berdasarkan ketiga metode diatas yang paling mendekati nilai exact dari ekspetasi nya adalah Metode Gauss Quadrature dengan nilai E(X)=0.7101422 pada iterasi ke-6 karena memiliki nilai error kecil dan hanya memerlukan 6 iterasi.

SOAL 2

Alt text

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-14

Metode Monte Carlo

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)
}

Perbandingan Metode untuk 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=50)
e1 = abs(exact_value -ht)

# Metode Simpson
hs  = simpson_n(fp,a=0,b=4,n=50)
e2 = abs(exact_value -hs)

# Metode Gauss Quadrature
gL  = gaussLegendre(n = 50,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(123)
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(50,50,50,m)
result  = c(ht,hs,hg,hmc)
error  = c(e1,e2,e3,e4)

library(DT)
Perbandingan = data.frame(Metode,n,result,error)
datatable(Perbandingan)

Berdasarkan perbandingan diatas dapat dilihat bahwa Metode Four-Gaus Quadrature merupakan metode terbaik dari ke-4 metode tersebut karena memiliki nilai error terkecil.

Referensi

Dito, G.A. (November 02, 2021). Metode Integral Numerik dengan R. Retrieved from https://gerrydito.github.io/Metode-Integral-Numerik-dengan-R/