runge-kutta method

Metode Runge-Kutta merupakan kelas metode yang secara bijak memanfaatkan informasi kemiringan kurva pada lebih dari satu titik untuk mengekstrapolasi solusi pada langkah selanjutnya [1]. Metode ini mengintegrasi secara numerik persamaan diferensial biasa menggunakan langkah coba-coba pada titik tengah suatu rentang agar suku kesalahan pada orde lebih rendah saling meniadakan [2]. Anggota keluarga metode Runge-Kutta yang paling dikenal adalah RK4, metode klasik Runge-Kutta atau cukup metode Runge-Kutta [3]. Telah terdapat contoh penerapan metode ini dalam berbagai bahasa pemrograman [4]. Metode ini pun dapat dijalankan beberapa processor secara paralel dalam menyelesaikan persamaan diferensial [5].

ode 1st order#

Suatu persamaan diferensial biasa (PDB) orde satu dapat dituliskan dalam bentuk

$$\begin{array}{}\text{(1)}& \frac{dy}{dx}=f(x,y)\end{array}$$

dengan $y=y(x)$.

euler method#

Dengan metode Euler Persamaan $\text{(1)}$ memiliki

$$\begin{array}{}\text{(2)}& y_{n+1}=y_n+hf(x_n,y_n)+O(h^2),\end{array}$$

sebagai solusi numeriknya [6], yang merupakan suatu permasalahan nilai awal, dengan $h$ adalah rentang antar langkah

$$\begin{array}{}\text{(3)}& x_{n+1}=x_n+h.\end{array}$$

Suku $O(h^2)$ pada Persamaan $\text{(2)}$ merupakan order akuras [7], yang dalam hal ini sebanding dengan $h^2$, dikenal sebagai simbol Landau [8].

rk2#

Formula untuk orde dua metode Runge-Kutta memiliki bentuk

$$\begin{array}{}\text{(4)}& k_1=hf(x_n,y_n),\end{array}$$

$$\begin{array}{}\text{(5)}& k_2=hf(x_{n+\frac{1}{2}h},y_{n+\frac{1}{2}{k}_1}),\end{array}$$

$$\begin{array}{}\text{(6)}& y_{n+1}=y_n+k_2+O(h^3),\end{array}$$

dan perlu ditambahkan pula Persamaan $\text{(3)}$.

rk4#

Formula untuk orde dua metode Runge-Kutta memiliki bentuk

$$\begin{array}{}\text{(7)}& k_1=hf(x_n,y_n),\end{array}$$

$$\begin{array}{}\text{(8)}& k_2=hf(x_{n+\frac{1}{2}h},y_{n+\frac{1}{2}{k}_1}),\end{array}$$

$$\begin{array}{}\text{(9)}& k_3=hf(x_{n+\frac{1}{2}h},y_{n+\frac{1}{2}{k}_2}),\end{array}$$

$$\begin{array}{}\text{(10)}& k_4=hf(x_{n+h},y_{n+k_3}),\end{array}$$

$$\begin{array}{}\text{(11)}& y_{n+1}=y_n+\frac{1}{6}{k}_1+\frac{1}{3}{k}_2+\frac{1}{3}{k}_3+\frac{1}{6}{k}_4+O(h^5),\end{array}$$

dan perlu ditambahkan pula Persamaan $\text{(3)}$.

example 1#

Misalkan terdapat suatu PDB orde pertama berbentuk

$$\begin{array}{}\text{(12)}& \frac{dy}{dx}=-cy\end{array}$$

yang akan digunakan sebagai contoh. Di sini $f(x,y)=f(y)$ dan $c$ suatu koefisien.

analytical solution#

Persamaan $\text{(12)}$ memiliki

$$\begin{array}{}\text{(13)}& y=A{e}^{-cx}\end{array}$$

sebagai solusi analitiknya dengan $A$ suatu koefisien.

euler solution#

Dengan menerapkan Persaamaan $\text{(2)}$ pada Persamaan $\text{(12)}$ akan diperoleh

$$\begin{array}{}\text{(14)}& \begin{array}{rcl}{y}_{n+1}& =& y_n-hc{y}_n\\ & =& (1-hc)y_n\end{array}\end{array}$$

sebagai solusinya dan dibutuhkan syarat awal $x_0$ dan $y_0=y(x_n)$ pada $n=0$.

rk2 solution#

Dengan menerapkan Persamaan $\text{(4)}$ pada Persamaan $\text{(12)}$ akan diperoleh

$$\begin{array}{}\text{(15)}& k_1=-hc{y}_n.\end{array}$$

Penggunaan Persamaan $\text{(5)}$ pada Persamaan $\text{(12)}$ menghasilkan

$$\begin{array}{}\text{(16)}& \begin{array}{rcl}{k}_2& =& -hc(y_n-\frac{1}{2}hc{y}_n)\\ & =& (-hc+\frac{1}{2}{h}^2{c}^2)y_n.\end{array}\end{array}$$

Dan substitusi Persamaan $\text{(15)}$ dan $\text{(16)}$ ke Persamaan $\text{(17)}$ akan memberikan

$$\begin{array}{}\text{(17)}& \begin{array}{rcl}{y}_{n+1}& =& y_n+(-hc+\frac{1}{2}{h}^2{c}^2)y_n\\ & =& (1-hc+\frac{1}{2}{h}^2{c}^2)y_n\end{array}\end{array}$$

sebagai solusinya dan dibutuhkan syarat awal $x_0$ dan $y_0=y(x_n)$ pada $n=0$.

rk4 solution#

Dengan menerapkan Persamaan $\text{(7)}$ pada Persamaan $\text{(12)}$ akan diperoleh

$$\begin{array}{}\text{(18)}& k_1=-hc{y}_n.\end{array}$$

Penggunaan Persamaan $\text{(8)}$ pada Persamaan $\text{(12)}$ menghasilkan

$$\begin{array}{}\text{(19)}& \begin{array}{rcl}{k}_2& =& -hc(y_n-\frac{1}{2}hc{y}_n)\\ & =& (-hc+\frac{1}{2}{h}^2{c}^2)y_n.\end{array}\end{array}$$

Substitusi Penggunaan Persamaan $\text{(9)}$ pada Persamaan $\text{(12)}$ memberikan

$$\begin{array}{}\text{(20)}& \begin{array}{rcl}{k}_3& =& -hc[y_n+\frac{1}{2}(-hc+\frac{1}{2}{h}^2{c}^2)y_n]\\ & =& (-hc+\frac{1}{2}{h}^2{c}^2-\frac{1}{8}{h}^3{c}^3)y_n.\end{array}\end{array}$$

Selanjutnya penggunaan Penggunaan Persamaan $\text{(9)}$ pada Persamaan $\text{(12)}$ akan memperoleh

$$\begin{array}{}\text{(21)}& \begin{array}{rcl}{k}_4& =& -hc[y_n+(-hc+\frac{1}{2}{h}^2{c}^2-\frac{1}{8}{h}^3{c}^3)y_n]\\ & =& (-hc+h^2{c}^2-\frac{1}{2}{h}^3{c}^3+\frac{1}{8}{h}^4{c}^4)y_n.\end{array}\end{array}$$

Dan substitusi Persamaan $\text{(18)}$, $\text{(19)}$, $\text{(20)}$, dan $\text{(21)}$ ke Persamaan $\text{(11)}$ akan menghasilkan

$$\begin{array}{}\text{(22)}& \begin{array}{rcl}{y}_{n+1}& =& y_n+\frac{1}{6}(-hc)y_n+\frac{1}{3}(-hc+\frac{1}{2}{h}^2{c}^2)y_n\\ & & +\frac{1}{3}(-hc+\frac{1}{2}{h}^2{c}^2-\frac{1}{8}{h}^3{c}^3)y_n\\ & & +\frac{1}{6}(-hc+h^2{c}^2-\frac{1}{2}{h}^3{c}^3+\frac{1}{8}{h}^4{c}^4)y_n\\ & =& (1-hc+\frac{1}{2}{h}^2{c}^2-\frac{1}{8}{h}^3{c}^3+\frac{1}{48}{h}^4{c}^4)y_n\end{array}\end{array}$$

sebagai solusinya dan dibutuhkan syarat awal $x_0$ dan $y_0=y(x_n)$ pada $n=0$.

the solutions#

Permasalahan pada Persamaan$\text{(12)}$ memiliki solusi analitik yang diberikan oleh Persamaan $\text{(13)}$ dan solusi analitiknya oleh Persamaan $\text{(14)}$, $\text{(17)}$, dan $\text{(22)}$.

Tabel 1. Beberapa hal terkait dengan contoh pada Persamaan $\text{(12)}$.

Solusi-solusi numerik yang diberikan pada Tabel 1 tidak selalu perlu untuk dituliskan dalam bentuk eksplisitnya sehingg terlihat jelas $y_{n+1}=y_{n+1}(y_n)$ akan tetapi cukup dengan menerapkan Persamaan-persamaan $\text{(15)}$, $\text{(16)}$, $\text{(17)}$ untuk metode RK2 dan Persamaan-persamaan $\text{(18)}$, $\text{(19)}$, $\text{(20)}$, $\text{(21)}$, $\text{(22)}$ untuk metode RK4. Contoh formula eksplisit pada Tabel 1 diberikan naya sebagai ilustrasi keterkaitan antar metode numerik.

example 2#

Suatu sistem fisis yang mirip dengan Persamaan$\text{(12)}$ adalah benda berbentuk bola yang jatuh dalam fluida karena tarikan gaya gravitasi bumi. Arah positif $y$ diambil ke atas dengan sehingga percepatan gravitasi $\overrightarrow{g}=-g_0\hat{y}$. Terdapat tiga gaya yang bekerja pada bola.

1. Gaya gravitasi bumi $$\begin{array}{}\text{(23)}& {\overrightarrow{F}}_G=-\frac{1}{6}\pi D^3\rho g_0\hat{y},\end{array}$$ dengan $\rho$ massa jenis bola dan $D$ adalah diameter bola.
2. Gaya apung fluida $$\begin{array}{}\text{(24)}& {\overrightarrow{F}}_B=\frac{1}{6}\pi D^3{\rho }_0{g}_0\hat{y},\end{array}$$ dengan ${\rho }_0$ massa jenis fluida.
3. Gaya gesek fluida $$\begin{array}{}\text{(25)}& {\overrightarrow{F}}_D=-3\pi {\eta }_{0}D\overrightarrow{v},\end{array}$$ dengan ${\eta }_0$ adalah kekentalan fluida.

Hukum Newton II memberikan

$$\begin{array}{}\text{(26)}& \begin{array}{rcl}\sum \overrightarrow{F}& =& m\overrightarrow{a}\\ {\overrightarrow{F}}_G+{\overrightarrow{F}}_B+{\overrightarrow{F}}_D& =& {\displaystyle \frac{1}{6}\pi D^3\rho \frac{d\overrightarrow{v}}{dt}}\\ -\frac{1}{6}\pi D^3\rho g_0\hat{y}& & \\ +\frac{1}{6}\pi D^3{\rho }_0{g}_0\hat{y}& & \\ -3\pi {\eta }_{0}D\overrightarrow{v}& =& \\ -\frac{1}{6}\pi D^3(\rho -{\rho }_0)g_0\hat{y}& & \\ -3\pi {\eta }_{0}D\overrightarrow{v}& =& \\ {\displaystyle -\left(\frac{\rho -{\rho }_0}{\rho }\right)g_0\hat{y}-\frac{3\pi {\eta }_{0}D}{\frac{1}{6}\pi D^3\rho }\overrightarrow{v}}& =& {\displaystyle \frac{d\overrightarrow{v}}{dt}}\\ {\displaystyle -(1-\frac{{\rho }_0}{\rho })g_0\hat{y}-\frac{18{\eta }_0}{D^2\rho }\overrightarrow{v}}& =& \end{array}\end{array}$$

yang dapat dituliskan tanpa notasi vektornya menjadi

$$\begin{array}{}\text{(27)}& \begin{array}{rcl}{\displaystyle \frac{dv}{dt}}& =& {\displaystyle -\frac{18{\eta }_0}{D^2\rho }v-(1-\frac{{\rho }_0}{\rho })g_0}\\ & =& {\displaystyle -\left(\frac{18{\eta }_0}{D^2\rho }\right)[v+\frac{D^2\rho }{18{\eta }_0}(1-\frac{{\rho }_0}{\rho })g_0].}\end{array}\end{array}$$

Misalkan suatu konstanta

$$\begin{array}{}\text{(28)}& c=\frac{18{\eta }_0}{D^2\rho },\end{array}$$

sehingga Persamaan $\text{(27)}$ dapat disederhanakan penulisannya menjadi

$$\begin{array}{}\text{(29)}& \frac{dv}{dt}=-c[v+\frac{g_0}{c}(1-\frac{{\rho }_0}{\rho })].\end{array}$$

Dapat didefinisikan suatu besaran baru

$$\begin{array}{}\text{(30)}& u=v+\frac{g_0}{c}(1-\frac{{\rho }_0}{\rho }),\end{array}$$

yang akan membuat Persamaan $\text{(29)}$ menjadi

$$\begin{array}{}\text{(31)}& \frac{du}{dt}=-cu,\end{array}$$

yang bentuknya telah serupa dengan Persamaan $\text{(12)}$. Dengan demikian solusi pada Tabel 1 dapat diterapkan pada Persamaan $\text{(31)}$.

real parameters#

Minyak motor SAE 50 memiliki viskositas (dinamis) $0.630\mathrm{P}\mathrm{a}\cdot \mathrm{s}$ pada $20{}^{\circ }\mathrm{C}$ [9], massa jenis aluminium adalah $2710\mathrm{k}\mathrm{g}/{\mathrm{m}}^3$ [10], dan massa jenis air murni adalah $998.2\mathrm{k}\mathrm{g}/{\mathrm{m}}^3$ pada $20{}^{\circ }\mathrm{C}$ [11], serta standar gravitasi bernilai $9.80665\mathrm{m}/{\mathrm{s}}^2$ [12]. Bola-bola aluminum dapat berukuran ${0.5}^{"}$, $1^{"}$, ${1.25}^{"}$, ${1.5}^{"}$, $2^{"}$, ${2.5}^{"}$, $3$, dan ${3.5}^{"}$ [13], dengan konversi $1^{"}=2.54\mathrm{c}\mathrm{m}$ [14]. Data-data ini dapat digunakan untuk mencari solusi dari sistem fisis pada Persamaan $\text{(27)}$.

exer#

1. Tunjukkan bagaimana pengaruh dari nilai $h$ terhadap hasil yang diperoleh dari solusi-solusi yang diberikan pada Tabel 1.
2. Apakah satuan dari $c$ pada Persamaan $\text{(28)}$?

notes#

1. Nikos Drakos, R. Sureshkumar, Michael Zeltkevic (translator), “Runge-Kutta Methods”, 10.001: Numerical Solution of Ordinary Differential Equations, 15 Apr 1998, url https://web.mit.edu/10.001/Web/Course_Notes/Differential_Equations_Notes/node5.html [20220322].
2. Eric W. Weisstein, “Runge-Kutta Method”, from MathWorld–A Wolfram Web Resource, Wolfram Research, Inc., 2022, url https://mathworld.wolfram.com/Runge-KuttaMethod.html [20220322].
3. Wikipedia contributors, “Runge–Kutta methods”, Wikipedia, The Free Encyclopedia, 9 March 2022, 07:21 UTC, url https://en.wikipedia.org/w/index.php?oldid=1076075901 [20220322].
4. GeeksforGeeks, Sam007, vt_m, divyesh072019, piotrkopec amartyaghoshgfg, “Runge-Kutta 4th Order Method to Solve Differential Equation”, GeeksforGeeks, 21 Jan 2022, url https://www.geeksforgeeks.org/runge-kutta-4th-order-method-solve-differential-equation/ [20220322].
5. Iman Al Fajri, Hendra Mesra, Jeffry Kusuma, “Solving Ordinary Differential Equation Using Parallel Fourth Order Runge-Kutta Method With Three Processors”, Jurnal Matematika, Statistika & Komputasi, [J Mat Stat Kom], vol 17, no 3, p 349-356, May 2021, url https://doi.org/10.20956/j.v17i3.12490.
6. freeCodeCamp.org, “Euler’s Method Explained with Examples”, freeCodeCamp, 26 Jan 2022, url https://www.freecodecamp.org/news/eulers-method-explained-with-examples/ [20220322].
7. Wikipedia contributors, “Order of accuracy”, Wikipedia, The Free Encyclopedia, 30 January 2022, 12:46 UTC, url https://en.wikipedia.org/w/index.php?oldid=1068844657 [20220322].
8. Eric W. Weisstein, “Landau Symbols”, from MathWorld–A Wolfram Web Resource, Wolfram Research, Inc., 2022, url https://mathworld.wolfram.com/LandauSymbols.html [20220322].
9. “Motor Oils - Dynamic Viscosities”, The Engineering ToolBox, 2011, url https://www.engineeringtoolbox.com/dynamic-viscosity-motor-oils-d_1759.html [20220322].
10. “Density of Aluminium”, Thyssenkrupp Materials (UK) Ltd, 2022, url https://www.thyssenkrupp-materials.co.uk/density-of-aluminium.html [20220322].
11. “Pure Water Density Standard”, Merck KGaA, Darmstadt, Germany, 2022, url https://www.sigmaaldrich.com/ID/en/product/sial/denwat [20220322].
12. “standard acceleration of gravity”, The NIST Reference on Constant, Units, and Uncertainty, url https://physics.nist.gov/cgi-bin/cuu/Value?gn [20220322].
13. “Aluminum Solid Ball No Hole”, Liberty Brass Turning Company, 2019, url https://libertybrass.com/products/aluminum-solid-ball-no-hole [20220322].
14. “Convert inches to cm”, unitconverters.net, 2022, url https://www.unitconverters.net/length/inches-to-cm.htm [20220322].