Manual

Metode Trapezoidal

Metode Trapezoidal Gambar 1

Metode Trapezoidal Gambar 2

Metode Simpson

Metode Simpson

Metode Gaus Quadrature

Metode Four-point Gauss Quadrature Gambar 1

Metode Four-point Gauss Quadrature Gambar 2

Metode Four-point Gauss Quadrature Gambar 3

Metode Four-point Gauss Quadrature Gambar 4

Program di R

Nilai exact dari nilai harapan distibusi kumaraswamy

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

## [1] 0.7102273

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


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

f <-function(x) 15*x^5*(1-x^5)^2
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-05

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


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

f2 <-function(x) 15*x^5*(1-x^5)^2
while(err>tol){
  res_simp <- simpson_n(f2,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

library(pracma)

## Warning: package 'pracma' was built under R version 4.0.5

#Mendefinisikan fungsi f(x)
f3 <- function(x) {15*x^5*(1-x^5)^2}

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

#Koefisien dan gauss point
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

#menghitung nilai integral
I <- sum(Ci * f3(xi))
I

## [1] 0.6792999

#menghitung nilai integral dengan toleransi
library(pracma)
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 * f3(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

#nilai exact
exact_value=0.7102273 
  
#trapezoidal
n=23 
result=0.7101278
error=9.953088e-05
trapezoidal <- cbind(n,result,error)

#simpson
n=34
result=0.7103096
error=8.225045e-05
simpson <- cbind(n,result,error)

#four-point Gauss Quadrature
n=6
result=0.7101422
error=8.507726e-05
four_point_Quadrature <- cbind(n,result,error)

exact_value

## [1] 0.7102273

komparasi <- rbind(trapezoidal,simpson, four_point_Quadrature)
komparasi

##       n    result        error
## [1,] 23 0.7101278 9.953088e-05
## [2,] 34 0.7103096 8.225045e-05
## [3,]  6 0.7101422 8.507726e-05

CDF

CDF Distribusi Eksponensial

Exact Value

exact_value2 <- 1- exp(-2*4)
exact_value2

## [1] 0.9996645

Trapezoidal

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

f <-function(x) 2*exp(-2*x)
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=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

Simpson

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

f2 <-function(x) 2*exp(-2*x)
while(err>tol){
  res_simp <- simpson_n(f2,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

Gauss Quadrature

library(pracma)

#Mendefinisikan fungsi f(x)
f3 <- function(x){
  2*exp(-2*x)
}

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

#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 nilai integral
I <- sum(Ci * f3(xi))
I

## [1] 0.9996645

#menghitung integral dengan toleransi
exact_value=0.9996645
tol <-  0.0001
err <- 1
n = 100

while(err>tol){
  
  gL <- gaussLegendre(n = n,a = 0,4)

 Ci <- gL$w # koefisien
 xi <- gL$x # gauss point
  
  res_gl <- sum(Ci * f3(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=100, result=0.9996645, error=3.737209e-08

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

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

set.seed(123)
int_g <- mc_integral(g,a=0,b=4,m = 1000)
int_g

## [1] 1.001826

#nilai exact
exact_value2 = 0.9996645 
  
#trapezoidal
n=231
result=0.9997644
error=9.995e-05
trapezoidal <- cbind(n,result,error)

#simpson
n=100
result=0.9996648
error=2.646781e-07
simpson <- cbind(n,result,error)

#four-point Gauss Quadrature
n=100
result=0.9996645
error=3.737209e-08
four_point_Quadrature <- cbind(n,result,error)

#Monte Carlo
monte_carlo = 1.001826

exact_value2

## [1] 0.9996645

komparasi <- rbind(trapezoidal,simpson, four_point_Quadrature)
komparasi

##        n    result        error
## [1,] 231 0.9997644 9.995000e-05
## [2,] 100 0.9996648 2.646781e-07
## [3,] 100 0.9996645 3.737209e-08

monte_carlo

## [1] 1.001826