Aproksimasi Integral

Berapakah nilai dari integral tentu ∫ba 𝑓(𝑥)𝑑𝑥 ?
Untuk fungsi seperti polinom ataupun fungsi yang cukup mudah dicari integralnya, nilai di atas dapat dicari.
Namun bagaimana jika 𝑓(𝑥) adalah suatu fungsi implisit atau fungsi yang cukup sulit untuk dicari integralnya? Contoh 𝑓(𝑥)=𝑒^(−𝑥2 )
Ada beberapa metode aproksimasi yang dapat digunakan untuk memberikan nilai hampiran dari integral tentu tersebut.
Metode yang akan dibahas: Riemann, Newton-Cotes, Romberg.
Secara umum untuk metode aproksimasi integral menggunakan cara mempartisi selang [a,b] menjadi n sub-selang.
𝑎=𝑥0<𝑥1<…<𝑥(𝑛−1)<𝑥𝑛=𝑏
Selanjutnya adalah bergantung dengan masing-masing metode dalam menggunakan titik hasil dari setiap partisi.

Riemann

Diperoleh sebanyak n sub-selang [𝑥𝑖,𝑥(𝑖+1)], dengan 𝑖=0,1,2,…,𝑛−1
Dipilih satu titik 𝑐_𝑖 pada setiap subselang yang akan dikenakan pada fungsi.
Jumlahan Riemann didefinisikan sebagai:
Untuk Riemann Kiri, dipilih 𝑐𝑖=𝑥𝑖
Untuk Riemann Kanan, dipilih 𝑐𝑖=𝑥(𝑖+1)
Untuk Riemann Tengah, dipilih 𝑐𝑖=(𝑥𝑖+𝑥(𝑖+1))/2
Dipilih pembagian sub-selang sehingga setiap sub-selang memiliki jarak yang sama: 𝑥(𝑖+1)−𝑥𝑖=(𝑏−𝑎)/𝑛
f = fungsi
b = batas atas
a = batas atas
n = banyak partisi
𝑥𝑖= a + (𝑖-1)d
𝑥𝑖+1 = a + 𝑖d
(𝑥𝑖+𝑥(𝑖+1))/2 = (a + (𝑖-1)d) + (a + 𝑖d )/ 2

berikut kode riemann kiri

riemann_kiri = function(f,a,b,n){
  d = (b-a)/n
  sum=0
  for(i in 1:n){
    sum=sum+f(a+(i-1)*d)*d
  }
  return(sum)
}
f = function(x) 1/(x^2)
riemann_kiri(f,1,5,100)

## [1] 0.8194644

#perbandingan dengan integral
integrate(f,1,5)

## 0.8 with absolute error < 1.1e-07

berikut kode riemann kanan

riemann_kanan = function(f,a,b,n){
  d = (b-a)/n
  sum=0
  for(i in 1:n){
    sum=sum+f(a+i*d)*d
  }
  return(sum)
}
f = function(x) x^2
riemann_kanan(f,0,2,1000)

## [1] 2.670668

#perbandingan dengan integral
integrate(f,0,2)

## 2.666667 with absolute error < 3e-14

berikut kode riemann tengah

riemann_tengah = function(f,a,b,n){
  d = (b-a)/n
  sum=0
  for(i in 1:n){
    sum=sum+f(a+i*d)*d
  }
  return(sum)
}
f = function(x) 12/8*x^3 + 11/2*x^2 + x
riemann_tengah(f,2,7,1000)

## [1] 1532.93

#perbandingan dengan integral
integrate(f,2,7)

## 1531.042 with absolute error < 1.7e-11

Newton-Cotes

Aproksimasi integral menggunakan Newton-Cotes diperoleh dari aproksimasi fungsi f(x) menggunakan polinomial Lagrange. (Review kembali Metode Numerik)
f = fungsi
b = batas atas
a = batas atas

berikut kode trapezoid

trapezoid=function(f,a,b){
  return (((b-a)/2)*(f(a)+f(b)))
}
f=function(x) exp(x)
trapezoid(f,0,3)

## [1] 31.62831

integrate(f,0,3)

## 19.08554 with absolute error < 2.1e-13

berikut kode simpson 1/3

simpson13=function(f,a,b){
  return (((b-a)/6)*(f(a)+4*f((a+b)/2)+f(b)))
}
f=function(x) exp(x)
simpson13(f,0,3)

## [1] 19.50615

integrate(f,0,3)

## 19.08554 with absolute error < 2.1e-13

berikut kode simpson 3/8

simpson38=function(f,a,b){
  h = (b-a)/3
  return (((b-a)/8)*(f(a)+3*f(a+h)+3*f(b-h)+f(b)))
}
f=function(x) exp(x)
simpson38(f,0,3)

## [1] 19.27783

integrate(f,0,3)

## 19.08554 with absolute error < 2.1e-13

berikut kode boole

boole=function(f,a,b){
  h = (b-a)/4
  return (((b-a)/90)*(7*f(a)+32*f(a+h)+12*f((a+b)/2)+32*f(b-h)+7*f(b)))
}
f=function(x) exp(x)
boole(f,0,3)

## [1] 19.09102

integrate(f,0,3)

## 19.08554 with absolute error < 2.1e-13

Keempat metode Newton-Cotes dapat diaplikasikan ke dalam bentuk kelipatan (Komposit).
Diterapkan untuk sub-interval hasil dari partisi.
Contoh: Multiple Trapezoid
Untuk interval yang dibagi menjadi n sub-interval, maka jumlahannya adalah
Dijabarkan menjadi

berikut kode multiple trapezoid

MulTrap=function(f,a,b,n){
  d=(b-a)/n
  sum=f(a)
  for(i in 1:n-1){
    sum=sum+2*f(a+i*d)
  }
  sum=sum+f(a+n*d)
  return (d*sum/2)
}
f=function(x) exp(x)
MulTrap(f,0,3,1000)

## [1] 19.08855

integrate(f,0,3)

## 19.08554 with absolute error < 2.1e-13

Untuk metode Simpson-1/3 juga dapat diaplikasikan dalam bentuk kelipatan
Syaratnya, jumlah partisi haruslah merupakan kelipatan 2.
Jumlahkan bentuk Simpson-1/3 dari indeks genap

berikut kode multiple simpson 1/3

MulSimp1=function(f,a,b,n){
  d=(b-a)/n
  sum=f(a)
  for(i in 1:n-1){
    if(i %% 2 == 1){
      sum=sum+4*f(a+i*d)
    } else {
      sum=sum+2*f(a+i*d)
    }
  }
  sum=sum+f(a+n*d)
  return (d*sum/3)
}
f=function(x) exp(x)
MulSimp1(f,0,3,1000)

## [1] 19.08754

integrate(f,0,3)

## 19.08554 with absolute error < 2.1e-13

Untuk metode Simpson-3/8 dan Boole juga dapat diberlakukan ke dalam kelipatan.
Metode Simpson-3/8: Syaratnya partisinya harus kelipatan 3
Nilai yang digunakan dalam Simpson 3/8 adalah 𝑥(𝑖+3)−𝑥𝑖=3ℎ
Sementara untuk Metode Boole: Syarat partisinya harus kelipatan 4
Nilai yang digunakan dalam Boole adalah 𝑥(𝑖+4)−𝑥𝑖=4ℎ

Romberg

Metode Romberg adalah improvisasi dari metode Newton-Cotes
Gabungan antara Multiple Trapezoid dan Ekstrapolasi Richardson.
Notasi dari Romberg adalah 𝑅_(𝑖,𝑗)
Nilai dari 𝑅_(𝑖,1) diperoleh dari integrasi menggunakan metode Multiple Trapezoid dengan banyaknya partisi adalah 2^(𝑖−1)
Pada halaman sebelumnya sudah dibuat fungsi MulTrap, jalankan fungsi untuk setiap iterasinya yaitu MulTrap(f,a,b,2^i-1)
Kemudian, untuk iterasi pada indeks j, digunakan formula sebagai berikut:

romberg=function(f,a,b,n){
  R=matrix(NA,nrow=n,ncol=n)
  R[1,1]=MulTrap(f,a,b,1)
  for (i in 2:n){
    R[i,1]=MulTrap(f,a,b,2^(i-1))
    for (j in 2:i){
      R[i,j]=(4^(j-1)*R[i,j-1]-R[i-1,j-1])/(4^(j-1)-1)
    }
  }
  return(R)
}
f=function(x) exp(x)
romberg(f,0,3,5)

##          [,1]     [,2]     [,3]     [,4]    [,5]
## [1,] 34.62831       NA       NA       NA      NA
## [2,] 24.03669 20.50615       NA       NA      NA
## [3,] 20.72190 19.61696 19.55769       NA      NA
## [4,] 19.68367 19.33760 19.31897 19.31519      NA
## [5,] 19.32892 19.21067 19.20221 19.20035 19.1999

integrate(f,0,3)

## 19.08554 with absolute error < 2.1e-13

Latihan

Gunakan semua teknik yang sudah dipelajari untuk memberikan aproksimasi dari integral dibawah ini:

Latihan 1

f = function(x) exp(3*x)*sin(2*x)
riemann_kiri(f,0,pi/4,1000)

## [1] 2.584487

riemann_kanan(f,0,pi/4,1000)

## [1] 2.592773

riemann_tengah(f,0,pi/4,1000)

## [1] 2.592773

trapezoid(f,0,pi/4)

## [1] 4.14326

simpson13(f,0,pi/4)

## [1] 2.583696

simpson38(f,0,pi/4)

## [1] 2.585789

boole(f,0,pi/4)

## [1] 2.587968

MulTrap(f,0,pi/4,1000)

## [1] 2.58863

MulSimp1(f,0,pi/4,1000)

## [1] 2.588629

romberg(f,0,pi/4,5)

##          [,1]     [,2]     [,3]     [,4]     [,5]
## [1,] 4.143260       NA       NA       NA       NA
## [2,] 2.973587 2.583696       NA       NA       NA
## [3,] 2.684173 2.587701 2.587968       NA       NA
## [4,] 2.612463 2.588560 2.588617 2.588627       NA
## [5,] 2.594584 2.588624 2.588628 2.588629 2.588629

integrate(f,0,pi/4)

## 2.588629 with absolute error < 2.9e-14

Latihan 2

f = function(x) cos(x^2)
riemann_kiri(f,-1,1,1000)

## [1] 1.809047

riemann_kanan(f,-1,1,1000)

## [1] 1.809047

riemann_tengah(f,-1,1,1000)

## [1] 1.809047

trapezoid(f,-1,1)

## [1] 1.080605

simpson13(f,-1,1)

## [1] 1.693535

simpson38(f,-1,1)

## [1] 1.760901

boole(f,-1,1)

## [1] 1.812769

MulTrap(f,-1,1,1000)

## [1] 1.810128

MulSimp1(f,-1,1,1000)

## [1] 1.809769

romberg(f,-1,1,5)

##          [,1]     [,2]     [,3]     [,4]     [,5]
## [1,] 2.161209       NA       NA       NA       NA
## [2,] 2.080605 2.053736       NA       NA       NA
## [3,] 2.009215 1.985418 1.980864       NA       NA
## [4,] 1.926593 1.899053 1.893295 1.891905       NA
## [5,] 1.872203 1.854074 1.851075 1.850405 1.850242

integrate(f,-1,1)

## 1.809048 with absolute error < 2e-14

Latihan 3

f = function(x) x/(sqrt(x^2-4))
riemann_kiri(f,3,3.5,1000)

## [1] 0.6362441

riemann_kanan(f,3,3.5,1000)

## [1] 0.6361826

riemann_tengah(f,3,3.5,1000)

## [1] 0.6361826

trapezoid(f,3,3.5)

## [1] 0.6400461

simpson13(f,3,3.5)

## [1] 0.6362387

simpson38(f,3,3.5)

## [1] 0.6362248

boole(f,3,3.5)

## [1] 0.6362135

MulTrap(f,3,3.5,1000)

## [1] 0.6368842

MulSimp1(f,3,3.5,1000)

## [1] 0.6366606

romberg(f,3,3.5,5)

##           [,1]      [,2]      [,3]      [,4]     [,5]
## [1,] 1.3108665        NA        NA        NA       NA
## [2,] 0.9726008 0.8598455        NA        NA       NA
## [3,] 0.8041641 0.7480185 0.7405633        NA       NA
## [4,] 0.7201274 0.6921152 0.6883883 0.6875601       NA
## [5,] 0.6781550 0.6641642 0.6623008 0.6618867 0.661786

integrate(f,3,3.5)

## 0.6362133 with absolute error < 7.1e-15

Latihan 4

f = function(x) 1/(x*log(x))
riemann_kiri(f,exp(1),2*exp(1),1000)

## [1] 0.5269415

riemann_kanan(f,exp(1),2*exp(1),1000)

## [1] 0.5262368

riemann_tengah(f,exp(1),2*exp(1),1000)

## [1] 0.5262368

trapezoid(f,exp(1),2*exp(1))

## [1] 0.647654

simpson13(f,exp(1),2*exp(1))

## [1] 0.5321106

simpson38(f,exp(1),2*exp(1))

## [1] 0.5292544

boole(f,exp(1),2*exp(1))

## [1] 0.5268155

MulTrap(f,exp(1),2*exp(1),1000)

## [1] 0.5275892

MulSimp1(f,exp(1),2*exp(1),1000)

## [1] 0.5272557

romberg(f,exp(1),2*exp(1),5)

##           [,1]      [,2]      [,3]      [,4]      [,5]
## [1,] 1.6476540        NA        NA        NA        NA
## [2,] 1.0609964 0.8654439        NA        NA        NA
## [3,] 0.7856090 0.6938132 0.6823711        NA        NA
## [4,] 0.6538760 0.6099650 0.6043751 0.6031371        NA
## [5,] 0.5896629 0.5682585 0.5654781 0.5648607 0.5647106

integrate(f,exp(1),2*exp(1))

## 0.526589 with absolute error < 5.8e-15