Metode Trapezoidal

User defined function untuk metode trapezoidal

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

Metode Simpson

User defined function untuk 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)
}

Metode Gauss Quadrature

Pada metode ini digunakan nilai koefisien dan titik gauss yang berada pada domain [-1,1] dengan fungsi gausslendre yang tersedia pada package pracma. Sebelumnya, batas \([a,b]\) akan ditransformasi sedemikian sehingga menjadi \([-1,1]\)

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

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

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

\[ \int_{a}^{b}f(x)dx = \int_{-1}^{1}f \left( \frac{1}{2}[t(b-a)+a+b] \right)\frac{1}{2}(b-a)dt \]

gauss <- function (fx, n, a, b ){
  
#Legendre order = n
gL <- gaussLegendre(n = n, a = a,b)

#mendefinisikan koefisien dan gauss point

Ci <- gL$w # koefisien
xi <- gL$x # gauss point

#Menghitung integral 
I <- sum(Ci * fx(xi))

return (list(I=I,Ci=Ci,xi=xi))
}

Metode Monte Carlo

User defined function 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)
}

a. Secara Manual

Diketahui pdf dari distribusi Kusmarawamy

\[ f(x|a,b) = abx^{(a-1)} (1-x^{(a)})^ {(b-1)}; x \in (0,1) \] dan diketahui a = 5 , b = 3, maka didapatkan

\[ f(x|5,3) = 5 \cdot 3x^{5-1}(1-x^5)^{3-1} = 15x^{4} (1-x^{5})^ {2} \] Ditanyakan \(E(x)\) untuk fungsi distribusi kumaraswamy tersebut.

\[ E(X) = \int_{-\infty}^{\infty} x f(x)dx \] \[ E(X) = \int_{0}^{1} x \cdot 15x^{4} (1-x^{5})^ {2}dx \]

Sehingga, fungsi yang akan dicari didefinisikan sebagai :

\[ E(X)= \int_{0}^{1} 15x^{5} (1-x^{5})^ {2}dx \] untuk selanjutnya dalam proses perhitungan, \(f(x)\) yang digunakan adalah fungsi \(E(x)\)

b. Menggunakan R

Mendefinisikan fungsi E(x)

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

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

## [1] 0.7102273

Fun <- function (fx, a, b) {
  integrate(fx, a,b)
}

Fun(f, 0,1)

## 0.7102273 with absolute error < 7.9e-15

trapezoid(f,0,1,n = 4)

## [1] 0.6312869

library(pracma)
trapzfun(f,0,1, maxit = 4)

## $value
## [1] 0.7098069
## 
## $iter
## [1] 4
## 
## $rel.err
## [1] 0.005901199

exact_value= exact
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.07894039
## n=5, result=0.6738123, error=0.03641494
## n=6, result=0.6914928, error=0.01873451
## n=7, result=0.6997126, error=0.0105147
## n=8, result=0.7039057, error=0.006321576
## n=9, result=0.7062116, error=0.004015628
## n=10, result=0.7075597, error=0.002667566
## n=11, result=0.7083885, error=0.001838794
## n=12, result=0.7089199, error=0.001307412
## n=13, result=0.7092729, error=0.0009543875
## n=14, result=0.7095147, error=0.0007126217
## n=15, result=0.7096846, error=0.0005426467
## n=16, result=0.7098069, error=0.0004203769
## n=17, result=0.7098966, error=0.0003306343
## n=18, result=0.7099637, error=0.0002635799
## n=19, result=0.7100146, error=0.0002126742
## n=20, result=0.7100538, error=0.0001734721
## n=21, result=0.7100844, error=0.0001428915
## n=22, result=0.7101085, error=0.000118756
## n=23, result=0.7101278, error=9.950361e-05

simpson_n(f,0,1,n = 4)

## [1] 0.7683974

exact_value= exact
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.05817008
## n=5, result=0.7683974, error=0.05817008
## n=6, result=0.7497082, error=0.03948097
## n=7, result=0.7497082, error=0.03948097
## n=8, result=0.728112, error=0.0178847
## n=9, result=0.728112, error=0.0178847
## n=10, result=0.7188088, error=0.008581559
## n=11, result=0.7188088, error=0.008581559
## n=12, result=0.7147289, error=0.00450162
## n=13, result=0.7147289, error=0.00450162
## n=14, result=0.712782, error=0.002554739
## n=15, result=0.712782, error=0.002554739
## n=16, result=0.711774, error=0.001546689
## n=17, result=0.711774, error=0.001546689
## n=18, result=0.7112144, error=0.0009871029
## n=19, result=0.7112144, error=0.0009871029
## n=20, result=0.7108852, error=0.0006578924
## n=21, result=0.7108852, error=0.0006578924
## n=22, result=0.7106819, error=0.0004545901
## n=23, result=0.7106819, error=0.0004545901
## n=24, result=0.7105511, error=0.0003238091
## n=25, result=0.7105511, error=0.0003238091
## n=26, result=0.710464, error=0.0002367075
## n=27, result=0.710464, error=0.0002367075
## n=28, result=0.7104042, error=0.000176941
## n=29, result=0.7104042, error=0.000176941
## n=30, result=0.7103621, error=0.0001348572
## n=31, result=0.7103621, error=0.0001348572
## n=32, result=0.7103318, error=0.0001045471
## n=33, result=0.7103318, error=0.0001045471
## n=34, result=0.7103096, error=8.227772e-05

ft <- function(t){
  (15/2)*((1/2*(t+1))^5)*(1-(1/2*(t+1))^5)^2
}

domaingauss <- gauss(ft, 4, -1,1)
domaingauss

## $I
## [1] 0.6792999
## 
## $Ci
## [1] 0.3478548 0.6521452 0.6521452 0.3478548
## 
## $xi
## [1] -0.8611363 -0.3399810  0.3399810  0.8611363

gL <- gaussLegendre(n = 4,a = -1,1)

Ci <- gL$w # koefisien
Ci

## [1] 0.3478548 0.6521452 0.6521452 0.3478548

xi <- gL$x # gauss point
xi

## [1] -0.8611363 -0.3399810  0.3399810  0.8611363

I <- sum(Ci * ft(xi))
I

## [1] 0.6792999

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

while(err>tol){
  
  gL <- gaussLegendre(n = n,a = -1,1)
  
  res_gl <- sum(Ci * ft(xi))
 Ci <- gL$w # koefisien
 xi <- gL$x # gauss point
  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.03092735
## n=5, result=0.6792999, error=0.03092735
## n=6, result=0.707375, error=0.002852251
## n=7, result=0.7101422, error=8.504999e-05

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

domainasal <- gauss(fg, 4, 0,1)
domainasal

## $I
## [1] 0.6792999
## 
## $Ci
## [1] 0.1739274 0.3260726 0.3260726 0.1739274
## 
## $xi
## [1] 0.06943184 0.33000948 0.66999052 0.93056816

gL <- gaussLegendre(n = 4,a = 0,1)

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 * fg(xi))
I

## [1] 0.6792999

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

while(err>tol){
  
  gL <- gaussLegendre(n = n,a = 0,1)
  
  res_gl <- sum(Ci * fg(xi))
 Ci <- gL$w # koefisien
 xi <- gL$x # gauss point
  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.03092735
## n=5, result=0.6792999, error=0.03092735
## n=6, result=0.707375, error=0.002852251
## n=7, result=0.7101422, error=8.504999e-05

c. Metode Terbaik

Jika diamati dari ketiga metode yang digunakan, ketiganya mampu menghasilkan nilai harapan yang mendekati nilai sebenarnya (0.7102273) dengan error yang bervariasi. Dalam kasus ini Metode Gauss Quadrature merupakan metode terbaik yang menghasilkan error terkecil diantara metode lainnya.

Selain itu untuk batas error toleransi yang ditetapkan, metode gauss quadrature secara komputasi juga memiliki kemampuan yang lebih cepat dibandingkan metode trapezoidal maupun simpson. Pada batas tol 0.0001, metode gauss quadrature hanya memerlukan hingga n = 7.

Metode Trapezoidal

1. Menggunakan user defined function

Sehingga untuk F(x) diperoleh

trapezoid(fe,0,4,n=100)

## [1] 1.000198

2. Menggunakan fungsi trapzfun dari package pracma

library(pracma)
trapzfun(fe,0,4, maxit = 100)

## $value
## [1] 0.9996646
## 
## $iter
## [1] 14
## 
## $rel.err
## [1] 5.958465e-08

3. Menghitung hingga hasil integral mendekati hasil exact-nya dengan toleransi 0.0001

exact_value= 0.9996645
tol <-  0.0001
err <- 1
n = 100

while(err>tol){
  res_trap <- trapezoid(fe,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=100, result=1.000198, error=0.0005331349
## n=101, result=1.000187, error=0.0005226319
## n=102, result=1.000177, error=0.0005124362
## n=103, result=1.000167, error=0.000502536
## n=104, result=1.000157, error=0.00049292
## n=105, result=1.000148, error=0.0004835774
## n=106, result=1.000139, error=0.0004744979
## n=107, result=1.00013, error=0.0004656717
## n=108, result=1.000122, error=0.0004570896
## n=109, result=1.000113, error=0.0004487425
## n=110, result=1.000105, error=0.000440622
## n=111, result=1.000097, error=0.00043272
## n=112, result=1.00009, error=0.0004250287
## n=113, result=1.000082, error=0.0004175406
## n=114, result=1.000075, error=0.0004102487
## n=115, result=1.000068, error=0.0004031462
## n=116, result=1.000061, error=0.0003962266
## n=117, result=1.000054, error=0.0003894836
## n=118, result=1.000047, error=0.0003829112
## n=119, result=1.000041, error=0.0003765039
## n=120, result=1.000035, error=0.0003702561
## n=121, result=1.000029, error=0.0003641625
## n=122, result=1.000023, error=0.0003582181
## n=123, result=1.000017, error=0.0003524181
## n=124, result=1.000011, error=0.0003467579
## n=125, result=1.000006, error=0.0003412329
## n=126, result=1, error=0.000335839
## n=127, result=0.9999951, error=0.0003305719
## n=128, result=0.9999899, error=0.0003254278
## n=129, result=0.9999849, error=0.0003204029
## n=130, result=0.99998, error=0.0003154935
## n=131, result=0.9999752, error=0.000310696
## n=132, result=0.9999705, error=0.0003060072
## n=133, result=0.9999659, error=0.0003014237
## n=134, result=0.9999614, error=0.0002969424
## n=135, result=0.9999571, error=0.0002925604
## n=136, result=0.9999528, error=0.0002882747
## n=137, result=0.9999486, error=0.0002840824
## n=138, result=0.9999445, error=0.000279981
## n=139, result=0.9999405, error=0.0002759677
## n=140, result=0.9999365, error=0.0002720401
## n=141, result=0.9999327, error=0.0002681958
## n=142, result=0.9999289, error=0.0002644324
## n=143, result=0.9999252, error=0.0002607477
## n=144, result=0.9999216, error=0.0002571395
## n=145, result=0.9999181, error=0.0002536057
## n=146, result=0.9999146, error=0.0002501442
## n=147, result=0.9999113, error=0.0002467531
## n=148, result=0.9999079, error=0.0002434306
## n=149, result=0.9999047, error=0.0002401747
## n=150, result=0.9999015, error=0.0002369837
## n=151, result=0.9998984, error=0.0002338558
## n=152, result=0.9998953, error=0.0002307895
## n=153, result=0.9998923, error=0.0002277832
## n=154, result=0.9998893, error=0.0002248352
## n=155, result=0.9998864, error=0.000221944
## n=156, result=0.9998836, error=0.0002191083
## n=157, result=0.9998808, error=0.0002163266
## n=158, result=0.9998781, error=0.0002135976
## n=159, result=0.9998754, error=0.0002109198
## n=160, result=0.9998728, error=0.0002082921
## n=161, result=0.9998702, error=0.0002057133
## n=162, result=0.9998677, error=0.000203182
## n=163, result=0.9998652, error=0.0002006972
## n=164, result=0.9998628, error=0.0001982577
## n=165, result=0.9998604, error=0.0001958624
## n=166, result=0.999858, error=0.0001935102
## n=167, result=0.9998557, error=0.0001912002
## n=168, result=0.9998534, error=0.0001889313
## n=169, result=0.9998512, error=0.0001867026
## n=170, result=0.999849, error=0.0001845131
## n=171, result=0.9998469, error=0.0001823618
## n=172, result=0.9998447, error=0.000180248
## n=173, result=0.9998427, error=0.0001781708
## n=174, result=0.9998406, error=0.0001761292
## n=175, result=0.9998386, error=0.0001741225
## n=176, result=0.9998366, error=0.00017215
## n=177, result=0.9998347, error=0.0001702108
## n=178, result=0.9998328, error=0.0001683042
## n=179, result=0.9998309, error=0.0001664294
## n=180, result=0.9998291, error=0.0001645858
## n=181, result=0.9998273, error=0.0001627727
## n=182, result=0.9998255, error=0.0001609893
## n=183, result=0.9998237, error=0.0001592352
## n=184, result=0.999822, error=0.0001575095
## n=185, result=0.9998203, error=0.0001558117
## n=186, result=0.9998186, error=0.0001541413
## n=187, result=0.999817, error=0.0001524976
## n=188, result=0.9998154, error=0.00015088
## n=189, result=0.9998138, error=0.0001492881
## n=190, result=0.9998122, error=0.0001477212
## n=191, result=0.9998107, error=0.0001461789
## n=192, result=0.9998092, error=0.0001446606
## n=193, result=0.9998077, error=0.0001431658
## n=194, result=0.9998062, error=0.0001416941
## n=195, result=0.9998047, error=0.000140245
## n=196, result=0.9998033, error=0.000138818
## n=197, result=0.9998019, error=0.0001374127
## n=198, result=0.9998005, error=0.0001360286
## n=199, result=0.9997992, error=0.0001346653
## n=200, result=0.9997978, error=0.0001333224
## n=201, result=0.9997965, error=0.0001319995
## n=202, result=0.9997952, error=0.0001306962
## n=203, result=0.9997939, error=0.0001294122
## n=204, result=0.9997926, error=0.0001281469
## n=205, result=0.9997914, error=0.0001269002
## n=206, result=0.9997902, error=0.0001256715
## n=207, result=0.999789, error=0.0001244606
## n=208, result=0.9997878, error=0.0001232671
## n=209, result=0.9997866, error=0.0001220908
## n=210, result=0.9997854, error=0.0001209311
## n=211, result=0.9997843, error=0.000119788
## n=212, result=0.9997832, error=0.0001186609
## n=213, result=0.999782, error=0.0001175497
## n=214, result=0.999781, error=0.0001164541
## n=215, result=0.9997799, error=0.0001153737
## n=216, result=0.9997788, error=0.0001143083
## n=217, result=0.9997778, error=0.0001132575
## n=218, result=0.9997767, error=0.0001122212
## n=219, result=0.9997757, error=0.0001111991
## n=220, result=0.9997747, error=0.0001101908
## n=221, result=0.9997737, error=0.0001091962
## n=222, result=0.9997727, error=0.0001082151
## n=223, result=0.9997717, error=0.000107247
## n=224, result=0.9997708, error=0.000106292
## n=225, result=0.9997698, error=0.0001053496
## n=226, result=0.9997689, error=0.0001044197
## n=227, result=0.999768, error=0.0001035021
## n=228, result=0.9997671, error=0.0001025965
## n=229, result=0.9997662, error=0.0001017028
## n=230, result=0.9997653, error=0.0001008207
## n=231, result=0.9997644, error=9.995e-05

Metode Simpson

1. Menggunakan user defined function

simpson_n(fe,0,4,n = 100)

## [1] 0.9996648

2. Menghitung hingga hasil integral mendekati hasil exact-nya dengan toleransi 0.0001

exact_value= 0.9996645
tol <-  0.0001
err <- 1
n = 100

while(err>tol){
  res_simp <- simpson_n(fe,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=100, result=0.9996648, error=2.646781e-07

Metode Gauss Quadrature

1. Menggunakan R dengan fungsi gaussLegendre dengan menggunakan domain \([-1,1]\)

Mengubah domain fungsi integral ke [-1, 1]

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

dan

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

\[dx=2dt\] sehingga dengan domain [-1,1] \(f(x)\) didefinisikan

\[ \int_{0}^{4} 2exp ^{-2 x } dx = \int_{-1}^{1} 2exp ^{-2 (2[t+1])} 2dt = \int_{-1}^{1} 4exp ^{(-4t-4)}dt=\int_{-1}^{1} 4exp ^{(-4(t+1))}dt \]

mendefinisikan fungsi hasil transformasi \(f(t)\)

ft <- function(t){
  4*exp(-4*(t+1))
}

Legendre order 100

gL <- gaussLegendre(n = 100, a = -1,1)

Koefisien dan titik gauss pada [1-1]

Ci <- gL$w # koefisien
Ci

##   [1] 0.0007346345 0.0017093927 0.0026839254 0.0036559612 0.0046244501
##   [6] 0.0055884280 0.0065469485 0.0074990733 0.0084438715 0.0093804197
##  [11] 0.0103078026 0.0112251140 0.0121314577 0.0130259479 0.0139077107
##  [16] 0.0147758845 0.0156296211 0.0164680862 0.0172904606 0.0180959407
##  [21] 0.0188837396 0.0196530875 0.0204032326 0.0211334421 0.0218430024
##  [26] 0.0225312203 0.0231974232 0.0238409603 0.0244612027 0.0250575445
##  [31] 0.0256294029 0.0261762192 0.0266974592 0.0271926134 0.0276611982
##  [36] 0.0281027557 0.0285168543 0.0289030896 0.0292610841 0.0295904881
##  [41] 0.0298909796 0.0301622651 0.0304040795 0.0306161866 0.0307983790
##  [46] 0.0309504789 0.0310723374 0.0311638357 0.0312248843 0.0312554235
##  [51] 0.0312554235 0.0312248843 0.0311638357 0.0310723374 0.0309504789
##  [56] 0.0307983790 0.0306161866 0.0304040795 0.0301622651 0.0298909796
##  [61] 0.0295904881 0.0292610841 0.0289030896 0.0285168543 0.0281027557
##  [66] 0.0276611982 0.0271926134 0.0266974592 0.0261762192 0.0256294029
##  [71] 0.0250575445 0.0244612027 0.0238409603 0.0231974232 0.0225312203
##  [76] 0.0218430024 0.0211334421 0.0204032326 0.0196530875 0.0188837396
##  [81] 0.0180959407 0.0172904606 0.0164680862 0.0156296211 0.0147758845
##  [86] 0.0139077107 0.0130259479 0.0121314577 0.0112251140 0.0103078026
##  [91] 0.0093804197 0.0084438715 0.0074990733 0.0065469485 0.0055884280
##  [96] 0.0046244501 0.0036559612 0.0026839254 0.0017093927 0.0007346345

xi <- gL$x # gauss point
xi

##   [1] -0.99971373 -0.99849195 -0.99629513 -0.99312494 -0.98898440 -0.98387754
##   [7] -0.97780936 -0.97078578 -0.96281365 -0.95390078 -0.94405587 -0.93328854
##  [13] -0.92160930 -0.90902957 -0.89556164 -0.88121868 -0.86601469 -0.84996453
##  [19] -0.83308388 -0.81538924 -0.79689789 -0.77762791 -0.75759812 -0.73682809
##  [25] -0.71533812 -0.69314920 -0.67028302 -0.64676191 -0.62260886 -0.59784747
##  [31] -0.57250193 -0.54659701 -0.52015802 -0.49321079 -0.46578165 -0.43789740
##  [37] -0.40958529 -0.38087298 -0.35178853 -0.32236034 -0.29261719 -0.26258812
##  [43] -0.23230248 -0.20178986 -0.17108008 -0.14020314 -0.10918920 -0.07806858
##  [49] -0.04687168 -0.01562898  0.01562898  0.04687168  0.07806858  0.10918920
##  [55]  0.14020314  0.17108008  0.20178986  0.23230248  0.26258812  0.29261719
##  [61]  0.32236034  0.35178853  0.38087298  0.40958529  0.43789740  0.46578165
##  [67]  0.49321079  0.52015802  0.54659701  0.57250193  0.59784747  0.62260886
##  [73]  0.64676191  0.67028302  0.69314920  0.71533812  0.73682809  0.75759812
##  [79]  0.77762791  0.79689789  0.81538924  0.83308388  0.84996453  0.86601469
##  [85]  0.88121868  0.89556164  0.90902957  0.92160930  0.93328854  0.94405587
##  [91]  0.95390078  0.96281365  0.97078578  0.97780936  0.98387754  0.98898440
##  [97]  0.99312494  0.99629513  0.99849195  0.99971373

Menghitung integral

I <- sum(Ci * ft(xi))
I

## [1] 0.9996645

2. Menggunakan R dengan fungsi gaussLegendre dengan menggunakan domain asal \([0,4]\)

fungsi dari domain asal [0,4]

fg <- function(x){
  2*exp(-2*x)
}

Legendre order 4 dengan domain \([0,4]\)

gL <- gaussLegendre(n = 100,a = 0,4)

mendefinisikan koefisien dan gauss point

Ci <- gL$w # koefisien
Ci

##   [1] 0.001469269 0.003418785 0.005367851 0.007311922 0.009248900 0.011176856
##   [7] 0.013093897 0.014998147 0.016887743 0.018760839 0.020615605 0.022450228
##  [13] 0.024262915 0.026051896 0.027815421 0.029551769 0.031259242 0.032936172
##  [19] 0.034580921 0.036191881 0.037767479 0.039306175 0.040806465 0.042266884
##  [25] 0.043686005 0.045062441 0.046394846 0.047681921 0.048922405 0.050115089
##  [31] 0.051258806 0.052352438 0.053394918 0.054385227 0.055322396 0.056205511
##  [37] 0.057033709 0.057806179 0.058522168 0.059180976 0.059781959 0.060324530
##  [43] 0.060808159 0.061232373 0.061596758 0.061900958 0.062144675 0.062327671
##  [49] 0.062449769 0.062510847 0.062510847 0.062449769 0.062327671 0.062144675
##  [55] 0.061900958 0.061596758 0.061232373 0.060808159 0.060324530 0.059781959
##  [61] 0.059180976 0.058522168 0.057806179 0.057033709 0.056205511 0.055322396
##  [67] 0.054385227 0.053394918 0.052352438 0.051258806 0.050115089 0.048922405
##  [73] 0.047681921 0.046394846 0.045062441 0.043686005 0.042266884 0.040806465
##  [79] 0.039306175 0.037767479 0.036191881 0.034580921 0.032936172 0.031259242
##  [85] 0.029551769 0.027815421 0.026051896 0.024262915 0.022450228 0.020615605
##  [91] 0.018760839 0.016887743 0.014998147 0.013093897 0.011176856 0.009248900
##  [97] 0.007311922 0.005367851 0.003418785 0.001469269

xi <- gL$x # gauss point
xi

##   [1] 0.0005725465 0.0030160987 0.0074097305 0.0137501259 0.0220312095
##   [6] 0.0322449186 0.0443812830 0.0584284485 0.0743726915 0.0921984341
##  [11] 0.1118882597 0.1334229299 0.1567814037 0.1819408580 0.2088767101
##  [16] 0.2375626412 0.2679706230 0.3000709442 0.3338322402 0.3692215233
##  [21] 0.4062042152 0.4447441807 0.4848037630 0.5263438204 0.5693237649
##  [26] 0.6137016013 0.6594339688 0.7064761830 0.7547822796 0.8043050595
##  [31] 0.8549961348 0.9068059759 0.9596839602 1.0135784216 1.0684367005
##  [36] 1.1242051957 1.1808294166 1.2382540368 1.2964229473 1.3552793122
##  [41] 1.4147656239 1.4748237593 1.5353950363 1.5964202718 1.6578398389
##  [46] 1.7195937255 1.7816215928 1.8438628344 1.9062566352 1.9687420312
##  [51] 2.0312579688 2.0937433648 2.1561371656 2.2183784072 2.2804062745
##  [56] 2.3421601611 2.4035797282 2.4646049637 2.5251762407 2.5852343761
##  [61] 2.6447206878 2.7035770527 2.7617459632 2.8191705834 2.8757948043
##  [66] 2.9315632995 2.9864215784 3.0403160398 3.0931940241 3.1450038652
##  [71] 3.1956949405 3.2452177204 3.2935238170 3.3405660312 3.3862983987
##  [76] 3.4306762351 3.4736561796 3.5151962370 3.5552558193 3.5937957848
##  [81] 3.6307784767 3.6661677598 3.6999290558 3.7320293770 3.7624373588
##  [86] 3.7911232899 3.8180591420 3.8432185963 3.8665770701 3.8881117403
##  [91] 3.9078015659 3.9256273085 3.9415715515 3.9556187170 3.9677550814
##  [96] 3.9779687905 3.9862498741 3.9925902695 3.9969839013 3.9994274535

Menghitung integral

I <- sum(Ci * fg(xi))
I

## [1] 0.9996645

Menghitung hingga hasil integral mendekati hasil exact-nya dengan toleransi \(0.0001\)

exact_value= 0.9996645
tol <-  0.0001
err <- 1
n = 100

while(err>tol){
  
  gL <- gaussLegendre(n = n,a = 0,4)
  
  res_gl <- sum(Ci * fe(xi))
 Ci <- gL$w # koefisien
 xi <- gL$x # gauss point
  err <- abs(res_gl-exact_value)
  
  cat("n=",n,", result=",res_gl,", error=",err,"\n",sep = "")
  
  n=n+1
  if(n==1000){
    break
  }
  
}

## n=100, result=0.9996645, error=3.737209e-08

Metode Monte Carlo

1. Menggunakan fungsi user defined function

set.seed(123)
mc_integral(fe,a=0,b=4, m=100)

## [1] 0.9330439

Metode Terbaik

Berdasarkan keempat metode yang dilakukan, metode gauss quadrature memiliki nilai error terkecil jika dibandingkan metode lainnya. Sehingga dalam hal ini, metode gauss quadrature merupakan metode terbaik dalam melakukan perhitungan integral numerik. Sebelumnya juga dilakukan percobaan pada n yang berbeda namun tidak ditampilkan, dari beberapa percobaan yang dilakukan, semakin besar n yang digunakan, akan memberikan presisi nilai integral yang semakin baik (error semakin kecil)