Integrasi Numerik

Metode Trapezoidal

Metode Trapezoidal yang digunakan untuk mendekati \(\int_{a}^{b}f(x)dx\) didefinisikan sebagai

\[ \begin{aligned} \int_{a}^{b}f(x)dx \approx T_{n} = \frac{h}{2} [f(x_{0})+2f(x_{1})+2f(x_{2})+\ldots+2f(x_{n-1})+f(x_{n})] \end{aligned} \] dimana \(h=\frac{b-a}{n}\) dan \(x_{i}=a+ih\) jika \(n \rightarrow \infty\), maka hasil dari \(\int_{a}^{b}f(x)dx\) sama dengan nilai integral definit. Berikut adalah 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

Metode Simpson yang digunakan untuk mendekati \(\int_{a}^{b}f(x)dx\) didefinisikan sebagai

\[ \int_{a}^{b}f(x)dx \approx S_{n} = \frac{h}{3} [f(x_{0})+4f(x_{1})+2f(x_{2})+4f(x_{3})+2f(x_{4})+\ldots+4f(x_{n-1})+f(x_{n})] \]

dimana \(h=\frac{b-a}{n}\) dan koefisien dari metode simpson mengikuti pola berikut:

\[\underbrace{1,4,2,4,2,\ldots,4,2,4,1}_{n+1 \space \text{titik}}\] Berikut adalah 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

Bentuk umum dari Gauss quadrature adalah

\[\int_{a}^{b}f(x)dx \approx I= \sum_{i=1}^{n}C_{i}f(x_{i})\]

dimana koefisien \(C_{i}\) merupakan bobot dan \(x_{i}\) merupakan titik Gauss yang berada dalam interval \([-1,1]\).

Nilai Bobot \(C_{i}\) dan titik Gauss \(x_{i}\) dapat diperoleh saat domain integralnya adalah \([-1,1]\) dan \(f(x)\) merupakan polynomial berderajat \(n\). Nilai-nilai tersebut secara lengkap bisa dilihat pada tabel berikut:

Nilai-nilai pada tabel ini juga bisa diakses dengan menggunakan fungsi gaussLegendre dari package pracma. Simbol w pada output dibawah ini merupakan bobot \(C_{i}\). Selain itu, nilai \(n\) pada fungsi gaussLegendre bisa lebih dari 6.

# untuk n=2
gaussLegendre(n=2,a=-1,b=1)

## $x
## [1] -0.5773503  0.5773503
## 
## $w
## [1] 1 1

# untuk n=3
gaussLegendre(n=3,a=-1,b=1)

## $x
## [1] -7.745967e-01  8.881784e-16  7.745967e-01
## 
## $w
## [1] 0.5555556 0.8888889 0.5555556

# untuk n=4
gaussLegendre(n=4,a=-1,b=1)

## $x
## [1] -0.8611363 -0.3399810  0.3399810  0.8611363
## 
## $w
## [1] 0.3478548 0.6521452 0.6521452 0.3478548

# untuk n=5
gaussLegendre(n=5,a=-1,b=1)

## $x
## [1] -9.061798e-01 -5.384693e-01  5.551115e-16  5.384693e-01  9.061798e-01
## 
## $w
## [1] 0.2369269 0.4786287 0.5688889 0.4786287 0.2369269

# untuk n=6
gaussLegendre(n=6,a=-1,b=1)

## $x
## [1] -0.9324695 -0.6612094 -0.2386192  0.2386192  0.6612094  0.9324695
## 
## $w
## [1] 0.1713245 0.3607616 0.4679139 0.4679139 0.3607616 0.1713245

# untuk n=7
gaussLegendre(n=7,a=-1,b=1)

## $x
## [1] -9.491079e-01 -7.415312e-01 -4.058452e-01  9.992007e-16  4.058452e-01
## [6]  7.415312e-01  9.491079e-01
## 
## $w
## [1] 0.1294850 0.2797054 0.3818301 0.4179592 0.3818301 0.2797054 0.1294850

Jika domain integral berupa \([a,b]\), maka perlu ada transformasi sedemikian sehingga domain integral nya menjadi \([-1,1]\)

\[\int_{a}^{b}f(x)dx \space \rightarrow \space \int_{-1}^{1}f(t)dt \] Transformasi ini dilakukan dengan mengubah variabel sebagai berikut

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

dan

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

Sehingga integralnya menjadi

\[ \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 \]

Metode Integral Monte Carlo

Algoritma Metode Monte Carlo

1. Tentukan fungsi \(f(x)\) dan selang integrasinya \([a,b]\)
2. Tentukan jumlah titik acak yang akan digunakan \(N\)
3. Lakukan produksi titik acak \(x\) dengan selang \([a,b]\) sejumlah \(N\)
4. Hitung \(f(x_i)\)
5. Hitung estimasi area

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

Package di R

Selain user-defined function di atas terdapat beberapa fungsi yang bisa langsung digunakan, yaitu

1. Metode Trapezoidal menggunakan fungsi trapzfun dari package pracma
2. Metode Simpson belum ada fungsi secara khusus
3. Metode Adaptive Quadrature menggunakan fungsi integrate
4. Metode Monte Carlo belum ada fungsi secara khusus.

Soal 1

Gunakan metode trapezoidal untuk mengintegralkan fungsi berikut untuk \(n=4\)

\[\int_{1}^{5} \frac{1}{x} dx \]

1. Hitunglah tanpa menggunakan software
2. Hitunglah menggunakan R dengan user-defined function
3. Hitunglah menggunakan R fungsi trapzfun dari package pracma
4. Hitunglah menggunakan R dengan user-defined function hingga hasil integral mendekati hasil exact-nya dengan toleransi \(0.0005\)

Jawab:

1. Hitunglah tanpa menggunakan software

Diketahui

\[f(x)=\frac{1}{x}, a=1,b=5\]

Panjang setiap sub interval adalah sebagai berikut:

\[h=\frac{b-a}{n}=\frac{5-1}{4}=\frac{4}{4}=1\] Karena \(n=4\) maka bentuk dari metode trapezoidal adalah

\[T_{4} = \frac{h}{2} [f(x_{0})+2f(x_{1})+2f(x_{2})+2f(x_{3})+f(x_{4})]\]

Menghitung fungsi \(f(x)\) pada setiap \(x_{i}\)

\[f(x_{0})=f(1)=\frac{1}{1}=1\] \[f(x_{1})=f(2)=\frac{1}{2}\]

\[f(x_{2})=f(3)=\frac{1}{3}\] \[f(x_{3})=f(4)=\frac{1}{4}\]

\[f(x_{4})=f(5)=\frac{1}{5}\]

karena \(h=1\) diperoleh

\[T_{4} = \frac{1}{2} [1+2 \cdot \frac{1}{2}+2 \cdot \frac{1}{3}+2 \cdot \frac{1}{4}+ \frac{1}{5}]=\frac{101}{60}=1.683333\] Sehingga

\[\int_{1}^{5} \frac{1}{x} dx \approx T_{4}=1.683333 \]

Jika kita bandingkan dengan hasil yang didapatkan secara exact adalah sebagai berikut:

\[\int_{1}^{5} \frac{1}{x} dx=\log{x}|_{1}^{5}=\log{5}-\log{1}=1.609\]

2. Hitunglah menggunakan R dengan user-defined function

Definisikan fungsi \(f(x)\)

f <- function(x){
  1/x
}

Menghitung integral menggunakan fungsi trapezoid

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

## [1] 1.683333

3. Hitunglah menggunakan R fungsi trapzfun dari package pracma

trapzfun(f, 1, 5, maxit =  4)

## $value
## [1] 1.614406
## 
## $iter
## [1] 4
## 
## $rel.err
## [1] 0.01456193

4. Hitunglah menggunakan R hingga hasil integral mendekati hasil exactnya dengan toleransi \(0.0005\)

exact_value=1.609
tol <-  0.0005
err <- 1
n = 4

while(err>tol){
  res_trap <- trapezoid(f,1,5,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=1.683333, error=0.07433333
## n=5, result=1.657907, error=0.04890706
## n=6, result=1.643601, error=0.03460084
## n=7, result=1.63478, error=0.02577994
## n=8, result=1.628968, error=0.01996825
## n=9, result=1.624942, error=0.01594152
## n=10, result=1.622039, error=0.01303901
## n=11, result=1.619879, error=0.01087908
## n=12, result=1.618229, error=0.009228993
## n=13, result=1.61694, error=0.007940381
## n=14, result=1.615915, error=0.006915077
## n=15, result=1.615086, error=0.006086061
## n=16, result=1.614406, error=0.005406324
## n=17, result=1.613842, error=0.004842114
## n=18, result=1.613369, error=0.004368695
## n=19, result=1.612968, error=0.003967606
## n=20, result=1.612625, error=0.003624844
## n=21, result=1.61233, error=0.003329638
## n=22, result=1.612074, error=0.003073587
## n=23, result=1.61185, error=0.002850067
## n=24, result=1.611654, error=0.002653798
## n=25, result=1.611481, error=0.002480525
## n=26, result=1.611327, error=0.002326794
## n=27, result=1.61119, error=0.002189776
## n=28, result=1.611067, error=0.002067133
## n=29, result=1.610957, error=0.001956925
## n=30, result=1.610858, error=0.001857527
## n=31, result=1.610768, error=0.00176757
## n=32, result=1.610686, error=0.001685896
## n=33, result=1.610612, error=0.001611519
## n=34, result=1.610544, error=0.001543595
## n=35, result=1.610481, error=0.0014814
## n=36, result=1.610424, error=0.001424306
## n=37, result=1.610372, error=0.001371771
## n=38, result=1.610323, error=0.001323323
## n=39, result=1.610279, error=0.001278548
## n=40, result=1.610237, error=0.001237084
## n=41, result=1.610199, error=0.001198614
## n=42, result=1.610163, error=0.001162854
## n=43, result=1.61013, error=0.001129558
## n=44, result=1.610099, error=0.001098503
## n=45, result=1.610069, error=0.001069494
## n=46, result=1.610042, error=0.001042353
## n=47, result=1.610017, error=0.001016925
## n=48, result=1.609993, error=0.0009930681
## n=49, result=1.609971, error=0.0009706553
## n=50, result=1.60995, error=0.0009495727
## n=51, result=1.60993, error=0.0009297169
## n=52, result=1.609911, error=0.0009109947
## n=53, result=1.609893, error=0.0008933214
## n=54, result=1.609877, error=0.0008766201
## n=55, result=1.609861, error=0.0008608208
## n=56, result=1.609846, error=0.0008458596
## n=57, result=1.609832, error=0.0008316785
## n=58, result=1.609818, error=0.0008182241
## n=59, result=1.609805, error=0.0008054475
## n=60, result=1.609793, error=0.000793304
## n=61, result=1.609782, error=0.0007817525
## n=62, result=1.609771, error=0.000770755
## n=63, result=1.60976, error=0.0007602769
## n=64, result=1.60975, error=0.0007502857
## n=65, result=1.609741, error=0.0007407519
## n=66, result=1.609732, error=0.0007316479
## n=67, result=1.609723, error=0.0007229484
## n=68, result=1.609715, error=0.0007146296
## n=69, result=1.609707, error=0.0007066697
## n=70, result=1.609699, error=0.0006990484
## n=71, result=1.609692, error=0.0006917466
## n=72, result=1.609685, error=0.0006847469
## n=73, result=1.609678, error=0.0006780327
## n=74, result=1.609672, error=0.0006715888
## n=75, result=1.609665, error=0.0006654008
## n=76, result=1.609659, error=0.0006594553
## n=77, result=1.609654, error=0.0006537399
## n=78, result=1.609648, error=0.0006482429
## n=79, result=1.609643, error=0.0006429532
## n=80, result=1.609638, error=0.0006378605
## n=81, result=1.609633, error=0.0006329552
## n=82, result=1.609628, error=0.0006282283
## n=83, result=1.609624, error=0.0006236711
## n=84, result=1.609619, error=0.0006192756
## n=85, result=1.609615, error=0.0006150343
## n=86, result=1.609611, error=0.0006109401
## n=87, result=1.609607, error=0.0006069861
## n=88, result=1.609603, error=0.0006031662
## n=89, result=1.609599, error=0.0005994743
## n=90, result=1.609596, error=0.0005959047
## n=91, result=1.609592, error=0.0005924521
## n=92, result=1.609589, error=0.0005891115
## n=93, result=1.609586, error=0.000585878
## n=94, result=1.609583, error=0.0005827471
## n=95, result=1.60958, error=0.0005797146
## n=96, result=1.609577, error=0.0005767763
## n=97, result=1.609574, error=0.0005739284
## n=98, result=1.609571, error=0.0005711672
## n=99, result=1.609568, error=0.0005684892
## n=100, result=1.609566, error=0.0005658912
## n=101, result=1.609563, error=0.0005633699
## n=102, result=1.609561, error=0.0005609224
## n=103, result=1.609559, error=0.0005585458
## n=104, result=1.609556, error=0.0005562374
## n=105, result=1.609554, error=0.0005539947
## n=106, result=1.609552, error=0.0005518151
## n=107, result=1.60955, error=0.0005496964
## n=108, result=1.609548, error=0.0005476362
## n=109, result=1.609546, error=0.0005456324
## n=110, result=1.609544, error=0.000543683
## n=111, result=1.609542, error=0.0005417861
## n=112, result=1.60954, error=0.0005399397
## n=113, result=1.609538, error=0.0005381422
## n=114, result=1.609536, error=0.0005363917
## n=115, result=1.609535, error=0.0005346867
## n=116, result=1.609533, error=0.0005330255
## n=117, result=1.609531, error=0.0005314068
## n=118, result=1.60953, error=0.0005298291
## n=119, result=1.609528, error=0.0005282909
## n=120, result=1.609527, error=0.0005267911
## n=121, result=1.609525, error=0.0005253282
## n=122, result=1.609524, error=0.0005239012
## n=123, result=1.609523, error=0.0005225089
## n=124, result=1.609521, error=0.00052115
## n=125, result=1.60952, error=0.0005198237
## n=126, result=1.609519, error=0.0005185288
## n=127, result=1.609517, error=0.0005172644
## n=128, result=1.609516, error=0.0005160295
## n=129, result=1.609515, error=0.0005148232
## n=130, result=1.609514, error=0.0005136446
## n=131, result=1.609512, error=0.0005124929
## n=132, result=1.609511, error=0.0005113673
## n=133, result=1.60951, error=0.000510267
## n=134, result=1.609509, error=0.0005091912
## n=135, result=1.609508, error=0.0005081392
## n=136, result=1.609507, error=0.0005071104
## n=137, result=1.609506, error=0.000506104
## n=138, result=1.609505, error=0.0005051193
## n=139, result=1.609504, error=0.0005041559
## n=140, result=1.609503, error=0.000503213
## n=141, result=1.609502, error=0.0005022901
## n=142, result=1.609501, error=0.0005013867
## n=143, result=1.609501, error=0.0005005021
## n=144, result=1.6095, error=0.0004996359

Soal 2

Gunakan metode simpson untuk mengintegralkan fungsi berikut untuk n=4, dengan

\[\int_{-2}^{2} 3^{x} dx \]

1. Hitunglah tanpa menggunakan software
2. Hitunglah menggunakan R
3. Hitunglah menggunakan R hingga hasil integral mendekati hasil exact-nya dengan toleransi \(0.0001\)

Jawab:

1. Hitunglah tanpa menggunakan software

Diketahui

\[f(x)=3^{x}, a=-2,b=2\]

Panjang setiap sub interval adalah sebagai berikut:

\[h=\frac{b-a}{n}=\frac{2-(-2)}{4}=\frac{4}{4}=1\]

Karena \(n=4\) maka bentuk dari metode Simpson adalah

\[S_{4} = \frac{h}{3} [f(x_{0})+4f(x_{1})+2f(x_{2})+4f(x_{3})+f(x_{4})]\]

Menghitung fungsi \(f(x)\) pada setiap \(x_{i}\)

\[f(x_{0})=f(1)=3^{-2}=\frac{1}{9}\] \[f(x_{1})=f(2)=3^{-1}=\frac{1}{3}\]

\[f(x_{2})=f(3)=3^{0}=1\] \[f(x_{3})=f(4)=3^{1}=3\]

\[f(x_{4})=f(5)=3^{2}=9\]

karena \(h=1\) diperoleh

\[S_{4} = \frac{1}{3} \left[\frac{1}{9}+4 \cdot \frac{1}{3}+2 \cdot 1+4 \cdot 3+9 \right]=\frac{220}{27}=8.148\] sehingga

\[\int_{-2}^{2} 3^{x} dx \approx S_{4}=8.148 \]

Jika kita bandingkan dengan hasil yang didapatkan secara exact adalah sebagai berikut:

\[\int_{-2}^{2} 3^{x} dx=\frac{3^{x}}{\log{3}} |_{-2}^{2}=\frac{3^{2}}{\log{3}}-\frac{3^{-2}}{\log{3}}=\frac{80}{9\log{3}}=8.091\]

2. Hitunglah menggunakan R

Definisikan fungsi \(f(x)\)

f2 <- function(x){
  3^{x}
}

Menghitung integral menggunakan metode simpson

simpson_n(f2,-2,2,n = 4)

## [1] 8.148148

3. Hitunglah menggunakan R hingga hasil integral mendekati hasil exact-nya dengan toleransi \(0.0001\)

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

while(err>tol){
  res_simp <- simpson_n(f2,-2,2,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=8.148148, error=0.05714815
## n=5, result=8.148148, error=0.05714815
## n=6, result=8.103168, error=0.01216778
## n=7, result=8.103168, error=0.01216778
## n=8, result=8.094965, error=0.003965356
## n=9, result=8.094965, error=0.003965356
## n=10, result=8.092654, error=0.001653864
## n=11, result=8.092654, error=0.001653864
## n=12, result=8.091811, error=0.0008110155
## n=13, result=8.091811, error=0.0008110155
## n=14, result=8.091447, error=0.0004466333
## n=15, result=8.091447, error=0.0004466333
## n=16, result=8.091269, error=0.0002688499
## n=17, result=8.091269, error=0.0002688499
## n=18, result=8.091174, error=0.0001739047
## n=19, result=8.091174, error=0.0001739047
## n=20, result=8.09112, error=0.0001195166
## n=21, result=8.09112, error=0.0001195166
## n=22, result=8.091087, error=8.656722e-05

Soal 3

Gunakan metode four-point gaussian quadratur untuk mengintegralkan fungsi berikut

\[\int_{0}^{3} \exp(-x^2) dx \]

1. Hitunglah tanpa menggunakan software

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

3. Hitunglah menggunakan R dengan fungsi gaussLegendre dengan menggunakan domain asal

4. Hitunglah menggunakan R hingga hasil integral mendekati hasil exact-nya dengan toleransi \(0.00001\)

5. Hitunglah tanpa menggunakan software

Diketahui

\[f(x)=\exp(-x^2), a=0,b=3\] Karena pada soal yang diinginkan adalah metode four-point gaussian quadratur, berarti nilai \(n=4\). Berdasarkan tabel gaussian quadratur nilai-nilai koefisien dan titik gauss adalah sebagai berikut:

gt::gt(data = bobot_gc[3,])

langkah 1

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]\] \[x=\frac{1}{2}[t(3-0)+3+0]\]

\[x=\frac{1}{2}[3t+3]=\frac{3}{2}[t+1]\]

dan

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

Sehingga integralnya menjadi

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

langkah 2

Kemudian, kita menggunakan nilai-nilai koefisien dan titik gauss saat \(n=4\)

\[ \int_{-1}^{1}\frac{3}{2}\exp\left(-\left[\frac{3}{2}[t+1]\right]^{2}\right)\approx I = C_{1}f(t_{1})+C_{2}f(t_{2})+C_{3}f(t_{3})+C_{4}f(t_{4}) \]

\[ I = C_{1}f(t_{1})+C_{2}f(t_{2})+C_{3}f(t_{3})+C_{4}f(t_{4}) \]

\[ I= 0.3478548f(-0.8611363)+0.6521452f(-0.3399810)+0.6521452f(0.3399810)+0.3478548f(0.8611363) \] Menghitung fungsi \(f(x)\) pada setiap \(x_{i}\)

\[f(-0.8611363)=\frac{3}{2}\exp\left(-\left[\frac{3}{2}[-0.8611363+1]\right]^{2}\right)=1.436311\]

\[f(-0.3399810)=\frac{3}{2}\exp\left(-\left[\frac{3}{2}[-0.3399810+1]\right]^{2}\right)=0.5628787\]

\[f(0.3399810)=\frac{3}{2}\exp\left(-\left[\frac{3}{2}[0.3399810+1]\right]^{2}\right)=0.02639659\]

\[f(0.8611363)=\frac{3}{2}\exp\left(-\left[\frac{3}{2}[0.8611363+1]\right]^{2}\right)=0.0006185\]

sehingga

\[ I = 0.347858 \cdot 1.436311+0.6521452 \cdot 0.5628787+0.6521452 \cdot 0.02639659+ 0.347858 \cdot 0.0006185 = 0.8841359 \] \[ \int_{0}^{3}\exp(-x^2)dx \approx I = 0.8841359 \] jika kita bandingkan dengan hasil yang didapatkan secara exact adalah sebagai berikut:

\[ \int_{0}^{3}\exp(-x^2)dx \] \[ = \frac{\sqrt\pi}{2} \int_{0}^{3}\frac{2\exp(-x^2)}{\sqrt\pi}dx \] bentuk \(\int_{0}^{z}\frac{2\exp(-x^2)}{\sqrt\pi}dx\) merupakan bentuk integral khusus yang bernama Gauss error function \(\text{erf}(z)\)

\[ =\frac{\sqrt\pi}{2} \text{erf}(3) = 0.886207 \]

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

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

ft <- function(t){
  (3/2)*exp((-(3/2*(t+1))^2))
}

Legendre order 4

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

mendefinisikan koefisien dan gauss point

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

Menghitung integral

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

## [1] 0.8841359

3. Hitunglah menggunakan R dengan fungsi gaussLegendre dengan menggunakan domain asal

Mendefinisikan fungsi \(f(x)\)

f3 <- function(x){
  exp(-(x^2))
}

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

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

Mendefinisikan koefisien dan gauss point

Ci <- gL$w # koefisien
Ci

## [1] 0.5217823 0.9782177 0.9782177 0.5217823

xi <- gL$x # gauss point
xi

## [1] 0.2082955 0.9900284 2.0099716 2.7917045

Menghitung integral

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

## [1] 0.8841359

4. Hitunglah menggunakan R hingga hasil integral mendekati hasil exact-nya dengan toleransi \(0.00001\)

exact_value=0.886207
tol <-  0.00001
err <- 1
n = 4

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

 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.8841359, error=0.00207107
## n=5, result=0.8865292, error=0.0003221738
## n=6, result=0.8861828, error=2.416802e-05
## n=7, result=0.8862082, error=1.22274e-06

Soal 4

Diketahui pdf (probability density function) dari distirbusi normal \(N(30,5)\) adalah sebagai berikut

\[ f(x|30,5)= \frac{1}{5\sqrt{2\pi}} \exp{\left[-\frac{1}{2} \left( \frac{x-30}{5} \right)^{2}\right]} \] Carilah cdf (cumulative density function) dari distribusi \(N(30,5)\) untuk \(x=32\) dengan menggunakan metode Monte Carlo dengan banyaknya sampling 1000.

Jawaban:

CDF dari suatu distribusi didefinisikan sebagai berikut

\[ F(x) = \int_{-\infty}^{x} f(t) dt \]

sehingga cdf dari distribusi \(N(30,5)\) untuk \(x=32\) adalah

\[ F(x=32|30,5) = \int_{-\infty}^{32} \frac{1}{5\sqrt{2\pi}} \exp{\left[-\frac{1}{2} \left( \frac{t-30}{5} \right)^{2}\right]} dt \]

karena batas bawah integral tersebut \(-\infty\), kita harus melakukan transformasi sedimikian sehingga batas bawah-nya terbatas di nilai tertentu. Menurut Chapra dan Canale (2015), jika batas atas dan bawah integral memiliki tanda berbeda (positive/negative) maka perlu dilakukan pembagian terlebih dahulu agar batas atas dan bawah integral memiliki tanda yang sama sebelum dilakukan transformasi

\[ F(x=32|30,5) = \int_{-\infty}^{-1} \frac{1}{5\sqrt{2\pi}} \exp{\left[-\frac{1}{2} \left( \frac{t-30}{5} \right)^{2}\right]} dt + \int_{-1}^{32} \frac{1}{5\sqrt{2\pi}} \exp{\left[-\frac{1}{2} \left( \frac{t-30}{5} \right)^{2}\right]} dt \]

Batas \(-1\) bisa diganti dengan bilangan negative apapun (dianjurkan untuk mendekati 0)

Kemudian integral

\[ \int_{-\infty}^{-1} \frac{1}{5\sqrt{2\pi}} \exp{\left[-\frac{1}{2} \left( \frac{t-30}{5} \right)^{2}\right]} dt \]

yang akan kita lakukan transformasi.

misalkan

\[ s = \frac{1}{t} \]

\[ ds =\frac{-1}{t^{2}} dt \]

\[ dt =-t^{2} ds \]

\[ dt =-\left(\frac{1}{s}\right)^{2} ds \] \[ dt =-\frac{1}{s^{2}} ds \]

sehingga domain integralnya pun berubah

\[t=-\infty \rightarrow s=\frac{1}{t}=\frac{1}{-\infty}=0 \] \[t=-1 \rightarrow s=\frac{1}{t}=\frac{1}{-1}=-1 \] jadi integral yang akan kita cari adalah

\[ \int_{0}^{-1} \frac{1}{5\sqrt{2\pi}} \exp{\left[-\frac{1}{2} \left( \frac{\frac{1}{s}-30}{5} \right)^{2}\right]} \left(-\frac{1}{s^{2}}\right) ds \]

batas atas dan bawah integral atas bisa dipertukarkan dengan mengalikan dengan \(-1\) agar memudahkan perhitungan

\[ \int_{-1}^{0} \frac{1}{5\sqrt{2\pi}} \exp{\left[-\frac{1}{2} \left( \frac{\frac{1}{s}-30}{5} \right)^{2}\right]} \left(\frac{1}{s^{2}}\right) ds \]

Selanjutnya kita akan menghitung cdf dengan

\[ F(x=32|30,5) = \int_{-1}^{0} \frac{1}{5\sqrt{2\pi}} \exp{\left[-\frac{1}{2} \left( \frac{\frac{1}{s}-30}{5} \right)^{2}\right]} \left(\frac{1}{s^{2}}\right) ds + \int_{-1}^{32} \frac{1}{5\sqrt{2\pi}} \exp{\left[-\frac{1}{2} \left( \frac{t-30}{5} \right)^{2}\right]} dt \]

Membuat fungsi di R berdasarkan hasil diatas

g1 <- function(s) {
  (1/(5*sqrt(2*pi)))*exp(-(1/2)*(((1/s)-30)/5)^(2))*(1/(s^2))
}

g2 <- function(x) {
  (1/(5*sqrt(2*pi)))*exp(-(1/2)*((x-30)/5)^(2))
}

set.seed(123)
int_g1 <- mc_integral(g1,a=-1,b=0,m = 1000)
int_g1

## [1] 2.841771e-10

set.seed(123)
int_g2 <- mc_integral(g2,a=-1,b=32,m = 1000)
int_g2

## [1] 0.6437677

int_g <- int_g1 + int_g2
int_g

## [1] 0.6437677

jadi nilai cdf dari \(N(30,5)\) untuk \(x=32\) adalah

\[ F(x=32|30,5) \approx 0.6437677 \]