Derivation of RK4

We first obtain $k_2$, $k_3$, and $k_4$ of the RK4 method based on $k_1=f(t,y(t))$ and the Taylor series expansion of the two-variable function $f(t,y(t))$:

$$k_2 = f(t+c_2h,\,y+ha_{21}k_1)=f(t+c_2h,\,y+hc_2f)$$

$$= f+hc_2(f_t+f_yf) +\frac{(hc_2)^2}{2}(f_{tt}+2ff_{ty}+f^2f_{yy}) +\frac{(hc_2)^3}{3!}(f_{ttt}+3ff_{tty}+3f^2f_{tyy}+f^3f_{yyy})+O(h^4)$$ (279)

$$k_3 = f(t+c_3h,\,y+h(a_{31}k_1+a_{32}k_2))$$

$$= f+h(c_3f_t+(a_{31}k_1+a_{32}k_2)f_y) +\frac{h^2}{2}(c_3^2f_{tt}+2c_3(a_{31}k_1+a_{32}k_2)f_{ty}+(a_{31}k_1+a_{32}k_2)^2f_{yy})$$

$$+\frac{h^3}{3!}(c_3^3f_{ttt}+3c_3^2(a_{31}k_1+a_{32}k_2)f_{tty} +3c_3(a_{31}k_1+a_{32}k_2)^2f_{tyy}+(a_{31}k_1+a_{32}k_2)^3f_{yyy})+O(h^4)$$

$$= f+hT_1+\frac{h^2}{2}T_2+\frac{h^3}{3!}T_3+O(h^4)$$ (280)

where $T_1$, $T_2$, and $T_3$ are respectively the coefficients of the first, second, and third order terms of $h$. To find these coefficients, we first consider

$$a_{31}k_1+a_{32}k_2 = a_{31}f+a_{32} \left(f+hc_2(f_t+f_yf)+\frac{(hc_2)^2}{2}(f_{tt}+2ff_{ty}+f^2f_{yy}) \right)+O(h^3)$$

$$= (a_{31}+a_{32})f+a_{32} \left(hc_2(f_t+f_yf)+\frac{(hc_2)^2}{2}(f_{tt}+2ff_{ty}+f^2f_{yy}) \right)+O(h^3)$$

$$= c_3f+a_{32}c_2(f_t+f_yf)h+a_{32}c_2^2(f_{tt}+2ff_{ty}+f^2f_{yy})\frac{h^2}{2}+O(h^3)$$

$$= c_3f+a_{32}c_2(f_t+f_yf)h+O(h^2)$$

$$= c_3f+O(h)$$ (281)

We then substitute the expressions of $a_{31}k_1+a_{32}k_2$ into $T_1,\;T_2$ and $T_3$ to get:

$$T_1 = c_3f_t+(a_{31}k_1+a_{32}k_2)$$

$$= c_3f_t+\left(c_3f+a_{32}c_2(f_t+f_yf)h+a_{32}c_2^2(f_{tt}+2ff_{ty}+f^2f_{yy})\frac{h^2}{2}\right)f_y +O(h^3)$$

$$= c_3(f_t+f_yf)+\left(a_{32}c_2(f_t+f_yf)h+a_{32}c_2^2(f_{tt}+2ff_{ty}+f^2f_{yy})\frac{h^2}{2}\right)f_y +O(h^3)$$

$$= c_3(f_t+f_yf)+a_{32}c_2(f_t+f_yf)f_yh +a_{32}c_2^2(f_{tt}+2ff_{ty}+f^2f_{yy})f_y\frac{h^2}{2}+O(h^3)$$

$$T_2 = c_3^2f_{tt}+2c_3(a_{31}k_1+a_{32}k_2)f_{ty}+(a_{31}k_1+a_{32}k_2)^2f_{yy}$$

$$= c_3^2f_{tt}+2c_3\left(c_3f+a_{32}c_2(f_t+f_yf)h\right)f_{ty} +(c_3f+a_{32}c_2(f_t+f_yf)h)^2f_{yy}+O(h^2)$$

$$= c_3^2f_{tt}+2c_3^2ff_{ty}+2a_{32}c_2c_3(f_t+f_yf)f_{ty}h +(c_3^2f^2+2c_3fa_{32}c_2(f_t+f_yf)h)f_{yy}+O(h^2)$$

$$= c_3^2(f_{tt}+2ff_{ty}+f^2f_{yy})+2a_{32}c_2c_3(f_t+f_yf)(f_{ty}+ff_{yy})h+O(h^2)$$

$$T_3 = c_3^3f_{ttt}+3c_3^2(a_{31}k_1+a_{32}k_2)f_{tty} +3c_3(a_{31}k_1+a_{32}k_2)^2f_{tyy}+(a_{31}k_1+a_{32}k_2)^3f_{yyy}$$

$$= c_3^3f_{ttt}+3c_3^2c_3f f_{tty}+3c_3c_3^2f^2f_{tyy}+c_3^3f^3f_{yyy}$$

$$= c_3^3(f_{ttt}+3ff_{tty}+3f^2f_{tyy}+f^3f_{yyy})+O(h)$$ (282)

Now we have

$$k_3 = f+hT_1+\frac{h^2}{2}T_2+\frac{h^3}{3!}T_3+O(h^4)$$

$$= f+h\left(c_3(f_t+f_yf)+a_{32}c_2(f_t+f_yf)f_yh +a_{32}c_2^2(f_{tt}+2ff_{ty}+f^2f_{yy})f_y\frac{h^2}{2}\right)$$

$$+\frac{h^2}{2}\left( c_3^2(f_{tt}+2ff_{ty}+f^2f_{yy})+2c_2c_3a_{32}(f_t+f_yf)(f_{ty}+ff_{yy})h\right)$$

$$+\frac{h^3}{3!}c_3^3(f_{ttt}+3ff_{tty}+3f^2f_{tyy}+f^3f_{yyy}) +O(h^4)$$

$$= f+c_3(f_t+f_yf)h+\left(a_{32}c_2(f_t+f_yf)f_y +\frac{1}{2}c_3^2(f_{tt}+2ff_{ty}+f^2f_{yy})\right)h^2$$

$$+\left(\frac{1}{2}a_{32}c_2^2(f_{tt}+2ff_{ty}+f^2f_{yy})f_y +c_2c_3a_{32} (f_t+f_yf)(f_{ty}+ff_{yy} )\right) h^3$$

$$+\frac{h^3}{3!}\left(c_3^3(f_{ttt}+3ff_{tty}+3f^2f_{tyy}+f^3f_{yyy})\right)+O(h^4)$$ (283)

$$k_4 = f(t+c_4h,\,y+h(a_{41}k_1+a_{42}k_2+a_{43}k_3))$$

$$= f+h(c_4f_t+(a_{41}k_1+a_{42}k_2+a_{43}k_3)f_y)$$

$$+ \frac{h^2}{2}(c_4^2f_{tt}+2c_4(a_{41}k_1+a_{42}k_2+a_{43}k_3)f_{ty}+(a_{41}k_1+a_{42}k_2+a_{43}k_3)^2f_{yy})$$

$$+ \frac{h^3}{3!}(c_4^3f_{ttt}+3c_4^2(a_{41}k_1+a_{42}k_2+a_{43}k_3)... ...a_{42}k_2+a_{43}k_3)^2f_{tyy}+(a_{41}k_1+a_{42}k_2+a_{43}k_3)^3f_{yyy}) +O(h^4)$$

$$= f+hT_1+\frac{h^2}{2}T_2+\frac{h^3}{3!}T_3+O(h^4)$$ (284)

where $T_1$, $T_2$, and $T_3$ are respectively the coefficients of the first, second, and third order terms of $h$. To find these coefficients, we first consider

$$a_{41}k_1+a_{42}k_2+a_{43}k_3$$

$$= a_{41}f+a_{42}\left(f+c_2(f_t+f_yf)h+c_2^2(f_{tt}+2ff_{ty}+f^2f_{... ..._t+f_yf)f_y+\frac{1}{2}c_3^2(f_{tt}+2ff_{ty}+f^2f_{yy})\right)h^2\right)+O(h^3)$$

$$= (a_{41}+a_{42}+a_{43})f+(a_{42}c_2+a_{43}c_3)(f_t+f_yf)h +(a_{42}... ...f_{tt}+2ff_{ty}+f^2f_{yy})\frac{h^2}{2} +a_{43}a_{32}c_2(f_t+f_yf)f_yh^2+O(h^3)$$

$$= c_4f+(a_{42}c_2+a_{43}c_3)(f_t+f_yf)h+O(h^2)$$

$$= c_4f+O(h)$$ (285)

We then substitute the expressions of $a_{41}k_1+a_{42}k_2+a_{43}k_3$ into $T_1,\;T_2$ and $T_3$ to get:

$$T_1 = c_4f_t+(a_{41}k_1+a_{42}k_2+a_{43}k_3)f_y$$

$$= c_4f_t+\left(c_4f+(a_{42}c_2+a_{43}c_3)(f_t+f_yf)h +(a_{42}c_2^2+... ..._{ty}+f^2f_{yy})\frac{h^2}{2} +a_{43}a_{32}c_2(f_t+f_yf)f_yh^2\right)f_y+O(h^3)$$

$$= c_4(f_t+f_yf)+(a_{42}c_2+a_{43}c_3)(f_t+f_yf)f_yh +\left(\frac{1}... ...2)(f_{tt}+2ff_{ty}+f^2f_{yy}) +a_{43}a_{32}c_2(f_t+f_yf)f_y\right)f_yh^2+O(h^3)$$

$$T_2 = c_4^2f_{tt}+2c_4(a_{41}k_1+a_{42}k_2+a_{43}k_3)f_{ty}+(a_{41}k_1+a_{42}k_2+a_{43}k_3)^2f_{yy}$$

$$= c_4^2f_{tt}+2c_4(c_4f+(a_{42}c_2+a_{43}c_3)(f_t+f_yf)h)f_{ty} +(c_4f+(a_{42}c_2+a_{43}c_3)(f_t+f_yf)h)^2f_{yy}+O(h^2)$$

$$= c_4^2f_{tt}+2c_4^2ff_{ty}+2c_4(a_{42}c_2+a_{43}c_3)(f_t+f_yf)h)f_{ty} +(c_4^2f^2+2c_4f(a_{42}c_2+a_{43}c_3)(f_t+f_yf)h)f_{yy}+O(h^2)$$

$$= c_4^2(f_{tt}+2ff_{ty}+f^2f_{yy})+2c_4(a_{42}c_2+a_{43}c_3)(f_t+f_yf)f_{ty}h +2c_4(a_{42}c_2+a_{43}c_3)(f_t+f_yf)ff_{yy}h+O(h^2)$$

$$= c_4^2(f_{tt}+2ff_{ty}+f^2f_{yy})+2c_4(a_{42}c_2+a_{43}c_3)(f_t+f_yf)(f_{ty} +f_{yy}f)h+O(h^2)$$

$$T_3 = c_4^3f_{ttt}+3c_4^2(a_{41}k_1+a_{42}k_2+a_{43}k_3)f_{tty} +3c_4(a_{41}^2k_1+a_{42}k_2+a_{43}k_3)^2f_{tyy}+(a_{41}k_1+a_{42}k_2+a_{43}k_3)^3f_{yyy}$$

$$= c_4^3f_{ttt}+3c_4^2c_4ff_{tty}+3c_4(c_4f)^2f_{tyy}+(c_4f)^3f_{yyy}+O(h)$$

$$= c_4^3(f_{ttt}+3ff_{tty}+3f^2f_{tyy}+f^3f_{yyy})+O(h)$$ (286)

Now we get

$$k_4 = f+hT_1+\frac{h^2}{2}T_2+\frac{h^3}{3!}T_3+O(h^4)$$

$$= f+h\left( c_4(f_t+f_yf)+(a_{42}c_2+a_{43}c_3)(f_t+f_yf)f_yh +\lef... ...)(f_{tt}+2ff_{ty}+f^2f_{yy}) +a_{43}a_{32}c_2(f_t+f_yf)f_y\right)f_yh^2 \right)$$

$$+\frac{h^2}{2} \left( c_4^2(f_{tt}+2ff_{ty}+f^2f_{yy})+2c_4(a_{42}c_2+a_{43}c_3)(f_t+f_yf)(f_{ty} +f_{yy}f)h \right)$$

$$+\frac{h^3}{3!}c_4^3(f_{ttt}+3ff_{tty}+3f^2f_{tyy}+f^3f_{yyy})+O(h^4)$$

$$= f+c_4(f_t+f_yf)h+(a_{42}c_2+a_{43}c_3)(f_t+f_yf)f_yh^2 +\left(\fr... ...43}c_3^2)(f_{tt}+2ff_{ty}+f^2f_{yy}) +a_{43}a_{32}c_2(f_t+f_yf)f_y\right)f_yh^3$$

$$+\frac{h^2}{2} c_4^2(f_{tt}+2ff_{ty}+f^2f_{yy})+c_4(a_{42}c_2+a_{43}c_3)(f_t+f_yf)(f_{ty} +f_{yy}f)h^3$$

$$+\frac{h^3}{3!}c_4^3(f_{ttt}+3ff_{tty}+3f^2f_{tyy}+f^3f_{yyy})+O(h^4)$$ (287)

We can now substitute the expressions for $k_1$, $k_2$, $k_3$, and $k_4$ into the RK4 iteration to get:

$$y(t+h) = y+h(b_1k_1+b_2k_2+b_3k_3+b_4k_4)$$

$$= y+hb_1f+hb_2\left[f+hc_2(f_t+f_yf)+\frac{(hc_2)^2}{2}(f_{tt}+2ff_... ...^2f_{yy}) +\frac{(hc_2)^3}{3!}(f_{ttt}+3ff_{tty}+3f^2f_{tyy}+f^3f_{yyy})\right]$$

$$+hb_3\left[f+c_3(f_t+f_yf)h+\left(a_{32}c_2(f_t+f_yf)f_y+\frac{1}... ...ht) h^3 +\frac{h^3}{3!} c_3^3(f_{ttt}+3ff_{tty}+3f^2f_{tyy}+f^3f_{yyy}) \right]$$

$$+hb_4\left[ f+c_4(f_t+f_yf)h+(a_{42}c_2+a_{43}c_3)(f_t+f_yf)f_yh^... ...yy}f)h^3 +\frac{h^3}{3!}c_4^3(f_{ttt}+3ff_{tty}+3f^2f_{tyy}+f^3f_{yyy}) \right]$$

$$+O(h^5)$$

$$= y+h(b_1+b_2+b_3+b_4)f+(b_2c_2+b_3c_3+b_4c_4)(f_t+f_yf)h^2 +\frac{h^3}{2}(b_2c_2^2+b_3c_3^2+b_4c_4^2)(f_{tt}+2ff_{ty}+f^2f_{yy})$$

$$+(b_3a_{32}c_2+b_4a_{42}c_2+b_4a_{43}c_3)(f_t+f_yf)f_yh^3 +\frac{... ...{2}(b_3a_{32}c_2^2+b_4a_{42}c_2^2+b_4a_{43}c_3^2)(f_{tt}+2ff_{ty}+f^2f_{yy})f_y$$

$$+\frac{h^4}{3!}(b_2c_2^3+b_3c_3^3+b_4c_4^3)(f_{ttt}+3ff_{tty}+3f^2f_{tyy}+f^3f_{yyy})$$

$$+(b_3c_3c_2a_{32}+b_4c_4c_2a_{42}+b_4c_4c_3a_{43})(f_t+f_yf)(f_{ty}+ff_{yy})h^4 +b_4a_{43}a_{32}c_2(f_t+f_yf)f_y^2h^4 +O(h^5)$$ (288)

On the other hand, the solution of the given differential equation $dy(t)/dt=f(t,y(t))$ can be Taylor expanded as shown below:

$$y(t+h) = y(t)+h\frac{dy}{dt}+\frac{h^2}{2!}\frac{d^2y}{dt^2} +\frac{h^3}{3!}\frac{d^3y}{dt^3}+\cdots+\frac{h^s}{4!}\frac{d^sy}{dt^s}+O(h^{s+1})$$

$$= y(t)+h f+\frac{h^2}{2!}\frac{df}{dt}+\frac{h^3}{3!}\frac{d^2f}{dt^2} +\frac{h^4}{4!}\frac{d^3f}{dt^3}+O(h^5)$$ (289)

Substituting the derivatives in Eq. (235) into this we get Eq. (241):

$$y(t+h) = y(t)+hf+O(h^2)$$

$$= y(t)+h f+\frac{h^2}{2!}(f_t+f_yf)+O(h^3)$$

$$= y(t)+h f+\frac{h^2}{2!}(f_t+f_yf) +\frac{h^3}{3!}(f_{tt}+2f_{ty}f +f_{yy}f^2+f_yf_t+f_y^2f)+O(h^4)$$

$$= y(t)+h f+\frac{h^2}{2!}(f_t+f_yf) +\frac{h^3}{3!}(f_{tt}+2f_{ty}f... ...h^3}{3!}f_y(f_t+f_yf) +\frac{h^4}{4!}(f_{ttt}+3f_{tty}f+3f_{tyy}f^2+f_{yyy}f^3)$$

$$+\frac{h^4}{4!}(2f_{ty}f_t+2f_{ty}f_yf+2f_{yy}f_tf+2f_{yy}f_yf^2+... ...{yy}f_tf+f_yf_{tt}+f_yf_{ty}f+2f_yf_{ty}f+2f_yf_{yy}f^2+f_y^2f_t+f_y^3f)+O(h^5)$$

$$= y(t)+h f+\frac{h^2}{2!}(f_t+f_yf) +\frac{h^3}{3!}(f_{tt}+2f_{ty}f... ...h^3}{3!}f_y(f_t+f_yf) +\frac{h^4}{4!}(f_{ttt}+3f_{tty}f+3f_{tyy}f^2+f_{yyy}f^3)$$

$$+\frac{h^4}{4!}f_y(f_{tt}+2f_{ty}f+f_{yy}f^2)+\frac{h^4}{4!}(3f_{ty}f_t+3f_{yy}f_tf+3f_{yy}f_yf^2+3f_yf_{ty}f+f_y^2f_t+f_y^3f)+O(h^5)$$

$$= y(t)+h f+\frac{h^2}{2!}(f_t+f_yf) +\frac{h^3}{3!}(f_{tt}+2f_{ty}f... ...h^3}{3!}f_y(f_t+f_yf) +\frac{h^4}{4!}(f_{ttt}+3f_{tty}f+3f_{tyy}f^2+f_{yyy}f^3)$$

$$+\frac{h^4}{24}f_y(f_{tt}+2f_{ty}f+f_{yy}f^2)+\frac{h^4}{8}(f_t+f_yf)(f_{ty}+f_{yy}f) +\frac{h^4}{24}(f_t+f_yf)f_y^2+O(h^5)$$ (290)

Finally, matching the coefficients in this Taylor expansion with those in the RK4 iteration obtained above, we get a set of necessary and sufficient conditions for the parameters of the RK4 method:

$$b_1+b_2+b_3+b_4=1$$

$$b_2c_2+b_3c_3+b_4c_4=1/2$$

$$b_2c_2^2+b_3c_3^2+b_4c_4^2=1/3$$

$$b_2c_2^3+b_3c_3^3+b_4c_4^3=1/4$$

$$b_3a_{32}c_2+b_4a_{42}c_2+b_4a_{43}c_3=1/6$$

$$b_3a_{32}c_3c_2+b_4a_{42}c_4c_2+b_4a_{43}c_4c_3=1/8$$

$$b_3a_{32}c_2^2+b_4a_{42}c_2^2+b_4a_{43}c_3^2=1/12$$

$$b_4a_{43}a_{32}c_2=1/24$$