INTEGRASI NUMERIK¶
Penyelesaian Integrasi dengan metode numerik ada beberapa macam, diantaranya adalah :
- Integrasi Reinmann
- Integrasi Trapezoida
- Integrasi Simpson
- Integrasi Romberg
Pada metode yang digunakan kali ini menggunakan Metode Romberg, Berikut penjelasan materi dari Metode tersebut.
Metode Romberg¶
Metode Romberg merupakan metode integrasi yang didasarkan pada perluasan ekstrapolasi Richardson yang dihasilkan dari aturan trapesium rekursif. Kelemahan dari metode ini adalah harus menggunakan jumlah interval yang besar guna mencapai akurasi yang diharapkan. Salah satu cara untuk meningkatkan akurasi adalah dengan membagi dua interval secara terus menerus sampai nilai integral yang dihitung dengan $2^k$ dan $2^{k+1}$ konvergen pada suatu nilai.
Metode ini sering digunakan untuk memperbaiki hasil aproksimasi oleh metode selisih terhingga. Metode ini dipakai untuk evalusasi numerik dari integrasl tertentu. Untuk dua interval bagian yang berbeda yang panjangnya $h_1$ dan $h_2$ akan diperoleh aproksimasi nilai-nilai $I_1$ dan $I_2$. Kemudian diperoleh kekeliruan $E_1$ dan $E_2$ $$ E_1 = -{1 \over 12} (b-a)h^2_1 y^{''}(x) \\ E_2 = -{1 \over 12} (b-a)h^2_2 y^{''}(x) $$ Karena suku $y^{''} (x=)$ adalah nilai terbesar dari $y^{''}(x)$ maka cukup beralasan untuk menganggap $y^{''}(x=) $ dan $y^{''}(x-)$ adalah sama. Karena $E_2-E_1 = I_1- I_2$, maka diperoleh : $$ E_2 = { h^2_2 \over h^2_2 h h^2_1} (I_2 - I_1) $$ Oleh karena itu aproksimasi baru diperoleh dengan bentuk : $$ I_3 = {I_1h^2_2 - I_2h^2_1 \over h^2_2 - h^2_1} $$ Karena menggunakan prinsip korektor, formula diatas akan memperbaiki nilai aproksimasi sebelumnya dan akan mendekati nilai yang sebenarnya. Dengan diperoleh persamaan baru yaitu : $$ I (h,{1 \over 2} h) = {1 \over 3}[4I({1 \over 2}h) - I(h)] \\ dengan \quad I(h) = I_1,I({1 \over 2}h) = I_2, dan \quad I(h,{1 \over 2}h) = I_3 $$ Metode integrasi Romberg didasarkan pada perluasan ekstrapolasi Richardson untuk memperoleh nilai integrasi yang semakin baik.
Proses integrasi Romberg
Kolom pertama pada tabel memuat hampiran integral tentu dengan menggunakan aturan trapesium rekursif. Kolom kedua merupakan hampiran integral yang sama dengan aturan Simpson rekursif (perbaikan pertama). Kolom ketiga merupakan hampiran integral yang sama dengan dengan aturan Boole rekursif (perbaikan kedua). Kolom keempat merupakan perbaikan ketiga. Demikian seterusnya
Metode Integrasi Romberg¶
Seperti halnya algoritma integrasi adaptif, integrasi Romberg adalah perluasan yang relatif mudah dari keluarga algoritma Newton-Cotes. Keduanya bekerja dengan menggunakan iterasi yang disempurnakan dari beberapa metode Newton-Cotes yang mendasarinya untuk memberikan perkiraan nilai integral yang lebih akurat. Tidak seperti proses komputasi fungsi riemann_adapint() , integrasi Romberg bukanlah pendekatan adaptif terhadap integrasi. Hal tersebut berarti metode Romberg tidak mengubah perilakunya sendiri berdasarkan perilaku fungsi yang akan diintegrasikan. Sebaliknya, kita mengeksploitasi perilaku fungsi trapesium pada batas untuk menghasilkan estimasi integral.
Pada Pemahaman integrasi Romberg, bisa kita mulai dengan implementasi rekursif dari aturan trapesium. Jik akita mulai dengna suatu fungsi , $T(f,m)$ dimana $T$ adalah fungsi trapesium, $f$ adalah fungsi yang akan diintegrasikan , dan $m$ adalah jumlah panel untuk diintegrasikan, maka : $$ S(f,m) = {4T(f,m) - T(f,m/2) \over 3} $$ dimana $S$ adalah aturan Simpson. Kemudian jika kita mendefinisikan $T=(f,0) = (b-a)(f(b)+f(a))=2,$ maka fungsi rekursif kita selesai, karena berdasarkan hubungan tersebut, fraksi yang diberikan dalam persamaan diatas juga merupakan perkiraan untuk integral.
Secara Umum : $$ I_{j,k} = {4^kI_{j,k-1} - I_{j-1,k-1} \over 4^k -1} $$ dimana $I_{0,0}$ adalah aturan trapsium satu panel dan $I_{j,0}$ adalah aturan trapesium dengan panel $2^j$. Dengan menggunakan fungsi - fungsi dasar ini $I_{j,l}$ dapat ditemukan secara iteratif sebagai matriks segitiga bawah dimana masing-masing nilai di kolom yang bukan paling kiri adalah fungsi dari nilai di sebelah kiri dan enti di atasnya.
Definisi rekursif ini muncul dari ekstrapolasi Richardson. Ketika diterapkan pada algoritma trapesium, yang konvergen menuju nilai sebenarnya dari integral sebagai $m$ (jumlah panel) meningkat, hubungan dalam Persamaan diatas muncul. Penting untuk dipahami bahwa pada batas ketika $k$ medekati tak terhingga , nilai $I_{j,k}$ adalah nilai sejati integral. Untuk nilai yang lebih kecil dari $k$, Integral Romberg hanya perkiraan, meskipun hasil yang diperoleh sangat bagus.
Pada Romberg Method
Proses pada Metode Romberg yang dimulai dari $R(0,0)$ berikut prosesnya
| R(0,0) | ||||
|---|---|---|---|---|
| R(1,0) | R(1,1) | |||
| R(2,0) | R(2,1) | R(2,2) | ||
| R(3,0) | R(3,1) | R(3,2) | R(3,3) | |
| ..... | ..... | .... | ..... | |
| R(n,0) | R(n,1) | R(n,2) | R(n,3) | R(n,m) |
Untuk $R(0,0)$ dapat kita definisikan sebagai berikut : $$ R(0,0) = {(b-a) \over 2} [f(a+f(b))], \\ R(n,0) = {1 \over 2}R(n-1,0) + h \sum^{2^{n+1}}_{k=1} f(a+(2k-1)h) \\ untuk \quad bisa \quad mencari \quad h, dengan: \\ h = {b-a \over 2^n} \\ R(n,m) = {1 \over 4^m - 1} [4^m X R(n,m-1) - R(n-1,m-1)] \quad syarat \quad n>=1,m>=1 \quad dan \quad m<n $$
Algoritma Penyelesaian Integrasi dengan Metode Romberg¶
-
Mendefinisikan fungsi integral $f(x)$
-
Menentukan batas-batas integral dengan nilai konstanta
-
Menentukan jumlah iterasi $(N)$
-
Menentukan nilai interval pada $[a,b]$ atau bisa dengan $[x_2, x_1]$ lalu menghitung h nya dengan cara : $h = x_2 -x_1$
-
Menghitung nilai integrasi pada kolom pertama dengan rumus : $$ R(1,1) = T_0 = {h \over 2}(f(x_1,y)+f(x_2,y)) $$
-
Menghitung nilai integrasi pada baris kedua sampai n pada kolom pertama dengan rumus : $$ R(n,1) = T_{k+1} = {T_k \over 2} + {h \over 2^{k+1}} {\sum^{2^k}{n=1}} {f{2r-1}} ; f(i) = f(x_1 + i {h \over 2^{k+1}}, n>=2) $$
-
Menghitung nilai integrasi pada kolom kedua sampai n dengan menggunakan rumus integrasi Romberg : $$ R(n,m) = {4^{m-1} R(n,m-1) - R(n-1,m-1) \over 4^{m-1}-1} $$
Implementasi Metode Romberg Dengan Python¶
Pada Integral yang diImplementasikan pada integral pada interval [0,1] of e^(-x^2), Berikut di Implementasikan dengan Python menggunakan import integrate
Listing Program
# import numpy and scipy.integrate import numpy as np from scipy import integrate # integral[0, 1] of e^(-x^2) interg = lambda x: np.exp(-x**2) # using scipy.integrate.romberg() romberg = integrate.romberg(interg, 0, 1, show = True) print(romberg)
Hasil Running (Output):
Romberg integration of <function vectorize1.<locals>.vfunc at 0x000001BF1D3A6F28> from [0, 1] Steps StepSize Results 1 1.000000 0.683940 2 0.500000 0.731370 0.747180 4 0.250000 0.742984 0.746855 0.746834 8 0.125000 0.745866 0.746826 0.746824 0.746824 16 0.062500 0.746585 0.746824 0.746824 0.746824 0.746824 32 0.031250 0.746764 0.746824 0.746824 0.746824 0.746824 0.746824 The final result is 0.7468241328122438 after 33 function evaluations. 0.7468241328122438
Pada Hasil Integrasi menggunakan fungsi def integrate dan tidak langsung menggunakan import integrate , juga menghitung tingkat error pada fungsi diatas yang telah dijalankan, berikut hasil codenya:
import numpy as np import sys def integrate(fn, a, b, steps=5, debug=False, exact=None): table = np.zeros((steps, steps), dtype=np.float64) pow_4 = 4 ** np.arange(steps, dtype=np.float64) - 1 # trapezoidal rule h = (b - a) table[0, 0] = h * (fn(a) + fn(b)) / 2 for j in range(1, steps): h /= 2 # extended trapezoidal rule table[j, 0] = table[j - 1, 0] / 2 table[j, 0] += h * np.sum(fn(a + i * h) for i in range(1, 2 ** j + 1, 2)) # richardson extrapolation for k in range(1, j + 1): table[j, k] = table[j, k - 1] + (table[j, k - 1] - table[j - 1, k - 1]) / pow_4[k] # debug if debug: print(table, file=sys.stderr) if exact is not None: errors = ['%.2e' % i for i in np.abs(table.diagonal() - exact)] print('abs. error:', errors, file=sys.stderr) return table[-1, -1] # integral[0, 1] of e^(-x^2) integrate(lambda x: np.exp(-x **2), 0, 1, debug=True, exact=0.746824132812427)
Dapat diketahui bahwa Nilai sebenarnya adalah 0.746824132812427, pada $e^{-x^2}$, Berikut Hasil Code yang telah dijalankan :
[[0.68393972058572 0. 0. 0. 0. ] [0.73137025182856 0.74718042890951 0. 0. 0. ] [0.74298409780038 0.74685537979099 0.74683370984975 0. 0. ] [0.7458656148457 0.74682612052747 0.7468241699099 0.74682401848228 0. ] [0.74658459678822 0.74682425743573 0.74682413322961 0.74682413264739 0.74682413309509]] abs. error: ['6.29e-02', '3.56e-04', '9.58e-06', '1.14e-07', '2.83e-10'] 0.7468241330950943
Contoh Lain menggunakan persamaan $f(x) = {1 \over (1+x)}$, berikut Hasil Code nya:
import numpy as np import sys def integrate(fn, a, b, steps=5, debug=False, exact=None): table = np.zeros((steps, steps), dtype=np.float64) pow_4 = 4 ** np.arange(steps, dtype=np.float64) - 1 # trapezoidal rule h = (b - a) table[0, 0] = h * (fn(a) + fn(b)) / 2 for j in range(1, steps): h /= 2 # extended trapezoidal rule table[j, 0] = table[j - 1, 0] / 2 table[j, 0] += h * np.sum(fn(a + i * h) for i in range(1, 2 ** j + 1, 2)) # richardson extrapolation for k in range(1, j + 1): table[j, k] = table[j, k - 1] + (table[j, k - 1] - table[j - 1, k - 1]) / pow_4[k] # debug if debug: print(table, file=sys.stderr) if exact is not None: errors = ['%.2e' % i for i in np.abs(table.diagonal() - exact)] print('abs. error:', errors, file=sys.stderr) return table[-1, -1] # integral[0, 1] of f(x) = 1/1+x print("Fungsi f(x) yang digunakan adalah") print("f(x) = 1/1+x") integrate(lambda x: 1/(1+x), 0, 1, debug=True, exact=2/3)
Hasil Outputnya :
Fungsi f(x) yang digunakan adalah f(x) = 1/1+x [[0.75 0. 0. 0. 0. ] [0.70833333333333 0.69444444444444 0. 0. 0. ] [0.69702380952381 0.69325396825397 0.6931746031746 0. 0. ] [0.69412185037185 0.69315453065453 0.69314790148123 0.69314747764483 0. ] [0.69339120220753 0.69314765281942 0.69314719429708 0.69314718307193 0.69314718191674]] abs. error: ['8.33e-02', '2.78e-02', '2.65e-02', '2.65e-02', '2.65e-02'] 0.693147181916745
Contoh Listing Program yang berbeda, Berikut Codenya
import numpy as np def trapezcomp(f, a, b, n): """ Composite trapezoidal function integration INPUTS: f: the function to integrate a: lower bound of integration b: upper bound n: number of panels to create between ``a`` and ``b`` """ # Initialization h = (b - a) / n x = a # Composite rule In = f(a) for k in range(1, n): x = x + h In += 2*f(x) return (In + f(b))*h*0.5 def romberg(f, a, b, p): """ Romberg integration INPUTS: f: the function to integrate a: lower bound of integration b: upper bound p: number of rows in the Romberg table """ I = np.zeros((p, p)) for k in range(0, p): # Composite trapezoidal rule for 2^k panels I[k, 0] = trapezcomp(f, a, b, 2**k) # Romberg recursive formula for j in range(0, k): I[k, j+1] = (4**(j+1) * I[k, j] - I[k-1, j]) / (4**(j+1) - 1) print("I[k, j+1] = (4**(j+1) * I[k, j] - I[k-1, j]) / (4**(j+1) - 1)=",I[k,0:k+1]) # display intermediate results return I if __name__ == '__main__': def func(x): return np.sin(x) p_rows = 4 I = romberg(func, 0, np.pi/2, p_rows) solution = I[p_rows-1, p_rows-1] print("Solusi yang Didapatkan = ",solution) # 1.00000000814
Hasil Output :
I[k, j+1] = (4**(j+1) * I[k, j] - I[k-1, j]) / (4**(j+1) - 1)= [0.78539816] I[k, j+1] = (4**(j+1) * I[k, j] - I[k-1, j]) / (4**(j+1) - 1)= [0.94805945 1.00227988] I[k, j+1] = (4**(j+1) * I[k, j] - I[k-1, j]) / (4**(j+1) - 1)= [0.9871158 1.00013458 0.99999157] I[k, j+1] = (4**(j+1) * I[k, j] - I[k-1, j]) / (4**(j+1) - 1)= [0.99678517 1.0000083 0.99999988 1.00000001] Solusi yang Didapatkan = 1.0000000081440203
Sekian TerimaKasih ;