## Fourth Order Runge Kutta Method Implemented on a Worksheet

The spreadsheet in Figure 10-2 illustrates the use of the RK method to simulate the first-order kinetic process A —> B, again using initial concentration [A]0 = 0.2000 and rate constant k = 5 x 10~3. The differential equation is, again, equation 10-4. This equation is of the simple form dy/dx = F(y), and thus only the yi terms of T\ to T4 need to be evaluated. The RK terms (note that T\ is the Euler method term) are shown in equations 10-12 through 10-15.

 t\ = -&[A], Ax (10-12) t2 = -£([A], + 7V2) Ax (10-13) ti - = -¿([A], + 7V2) Ax (10-14) t4 = -¿([A], + t3) ax (10-15)
 A | B C D E F G 1 Ruiiije-Kirtta simulation of firsl order kinetics 2 rate constant = 5.0E-03 « 3 time increment = 20

Figure 10-2. Simulation of first-order kinetics by the Runge-Kutta method, (folder 'Chapter 10 Examples', workbook 'ODE Examples', worksheet 'RK1')

Figure 10-2. Simulation of first-order kinetics by the Runge-Kutta method, (folder 'Chapter 10 Examples', workbook 'ODE Examples', worksheet 'RK1')

The RK equations in cells B7, C7, D7, E7 and F7, respectively, are (only part of the spreadsheet is shown; the formulas extend down to row 74):

If you use the names TAI, ..., TA4 you can use AutoFill to generate the column labels TAI, ..., TA4. These names are accepted by Excel, whereas T1 is not a valid name. As well, the nomenclature is expandable to systems requiring more than one set of Runge-Kutta terms (e.g., TBI, ..., TB4, etc.).

Compare the RK result in column F of Figure 10-2 with the analytical expression for the concentration, [A]/ = [A]0e ~k', in column G. After one half-life (row 13) the RK calculation differs from the analytical expression by only

0.00006%. (Compare this with the 3.6% error in the Euler method calculation at the same point.) Even after 10 half-lives (not shown), the RK error is only 0.0006%.

In essence, the fourth-order Runge-Kutta method performs four calculation steps for every time interval. The percent error after one half-life {t = 140) is only 6 x 10"5%. In contrast, in the solution by Euler's method, decreasing the time increment to 5 seconds to perform four times as many calculation steps still only reduces the error to 0.9% after 1 half-life.

If the spreadsheet is constructed as shown in Figure 10-2, you can't use a formula in which a name is assigned to the values of the calculated concentration in column F (the range \$F\$7:\$F\$74). This is because the formula in B7, for example, will use the concentration in F7; this is called an implicit intersection. An alternative arrangement that permits using a name for the concentration [A], is shown in Figure 10-3. Each row contains the concentration at the beginning and at the end of the time interval. The name C_t can now be assigned to the array of values in column B; the former formulas (now in cells \$C\$7:\$G\$74) contain C_t in place of F6 and cell B7 contains the formula =G6.

 A B C D E F G H 1 Runge-Kutta simulation of first or der kinetics 2 ■ -------- rate constant = 5.0E-03 (10 3 time increment = 20 (DX) 4 5 t Ct TAI -0.0190 TA3 TA4 RK exact 6 0 0.2000 -0,0200 -0.0191 -0.0181 0.1810 0.1810 7 20 0.1810 -0.01 81 -0.0172 -0.0172 -0.0164 0.1537 0.1637 8 40 0:1637 -0.0164 -0.0156 -0.0156 -0.0148 0.1482 0.1482 3 60 0.1482 -0,0148 -0.0141 -0.0141 -0.0134 0.1341 0.1341 10 80 0.1341 -0.01 34 -0.0127 -0.0128 -0.0121 0.1213 0.1213 11 100 0.1213 -0.0121 -0.0115 -0.0116 -0.0110 0.1038 0.1098 12 120 0.1098 -0.0110 -0.0104 -0.0105 -0.0039 0 0993 0.0993 13 140 0.0993 -0.0099 -0.0034 -0.0095 -0.0090 0,0899 0.0899 14 160 0.0899 -0.0090 -0.0085 -0.0086 -0.0081 0,0813 0.0813

Figure 10-3. Alternative spreadsheet layout for the Runge-Kutta method, (folder 'Chapter 10 Examples', workbook 'ODE Examples', worksheet 'RK2')

Figure 10-3. Alternative spreadsheet layout for the Runge-Kutta method, (folder 'Chapter 10 Examples', workbook 'ODE Examples', worksheet 'RK2')

The RK equations in cells C6, D6, E6, F6 and G6, respectively, are

=-k*(C_t+TA3)*DX =C_t+(TA1 +2*TA2+2*TA3+TA4)/6 and cell B7 contains the formula =G6.