## Fourth Order Runge Kutta Custom Function for Systems of Differential Equations

The simple Runge-Kutta custom function of Figure 10-4 was expanded so as to handle multiple differential equations, by using equations 10-21 through 1023. The VBA code is shown in Figure 10-9.

### The syntax of the custom function is

Runge3(x_var/aWe, y_variables, deriv_formulas, interval, index). The argument x_variable is a reference to the cell containing the independent variable, the argument y_variables is a reference to the range containing the values of the N dependent variables, and the argument deriv_formulas is a reference to the range containing the formulas of the N derivatives, in the same order as y_variables. For y_variables and deriv_formulas, the user can enter a range of cells or make a nonadjacent selection. The argument increment is the Ax used in the calculation. The optional argument index specifies the dependent variable to return; if omitted, the function returns the complete array of dependent variables. In this case the user must select a range of cells in a row, enter the formula and then press CONTROL+SHIFT+ENTER. Since the function always calculates the complete array, this can save calculation time if several dependent variables are being returned.

Option Explicit Option Base 1

Function Runge3(x_variable, y_variables, deriv_formulas, interval, Optional _ index)

'Runge-Kutta method to solve ordinary differential equations.

'Solves problems involving simultaneous first-order differential equations.

'x_variable is a reference to the independent variable x. 'y_variables is a reference to the dependent variables y(1)... y(N). 'deriv_formulas is a reference to the derivatives dy(i)/dx, in same order, 'interval is a reference to delta x

'index specifies the y(i) to be returned. If omitted, returns the array.

Dim FormulaText() As String, XAddr As String, YAddr() As String Dim J As Integer, N As integer

N = y_variables.Columns.Count

If N = 1 Then N = y_variables.Rows.Count

FormulaText(J) = Application.ConvertFormula(deriv_formulas(J).Formula, xlA1, xlA1, xlAbsolute) Next J

If IsMissing(index) Then

Runge3 = RK3(N, FormulaText, XAddr, YAddr, x_variable, y_variables, _ interval) (index) End If

End Function

I,. I,, I ,11. I.__, , , , , , , , , . . , I . . I . . I I

Dim X As Double, Y() As Double, term() As Double Dim J As Integer, K As Integer ReDim term(4, N), Y(N)

For J = 1 To N: Y(J) = y_variables(J).Value: Next J Call eval3(N, FormulaText, XAddr, YAddr, X, Y, H, K, term) K = 2: X = x_variable.Value + H / 2

For J = 1 To N: Y(J) = y_variables(J).Value + term(1, J) / 2: Next J

K = 3: X = x_variable.Value + H / 2 For J = 1 To N: Y(J) = y_variables(J) .Value + term(2, J) / 2: Next J Call eval3(N, FormulaText, XAddr, YAddr, X, Y, H, K, term) K = 4: X = x_variable.Value + H

For J = 1 To N: Y(J) = y_variables(J).Value + term(3, J): Next J Call eval3(N, FormulaText, XAddr, YAddr, X, Y, H, K, term)

Y(J) = y_variables(J).Value+(term(1, J)+2*term(2, J)+2*term(3, J)+term(4, J)) / 6 Next J RK3 = Y End Function

Sub eval3(N, FormulaText, XAddr, YAddr, X, Y, H, K, term) Dim I As Integer, J As Integer Dim T As String

T = FormulaText(J)

Call SubstitutelnString(T, YAddr(l), Y(l)) Next I

'+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Sub SubstitutelnString(T, Ref, Value)

'Replaces each instance of e.g., \$A\$2 in formula with number value, e.g., 0.20, then evaluates.

'Must do this replacement from end of formula to beginning.

'Modified 03/08/06 to handle possible un-intended replacement of e.g., \$A\$2 in

'Method: replace \$A\$2 with value & " "

'so that \$A\$22 becomes "0.20 2" and this formula evaluates to an error. Dim temp As String

Dim NReplacements As Integer, J As Integer Dim dummy As Double

'Substitute all instances of address with value

NReplacements = (Len(T) - Len(Application.Substitute(T, Ref, ""))) / Len(Ref) For J = NReplacements To 1 Step -1 temp = Application.Substitute^, Ref, Value & " ", J) On Error GoTo ErrorHandler dummy = Evaluate(temp) T = temp pt1: Next J Exit Sub

ErrorHandler:

'Trappable error number 13 (Type mismatch) is expected._

If Err.Number = 13 Then On Error GoTo 0 'Disable the error handler. Resume pt1 'and continue execution. Else

End 'Some other error, so quit completely End If End Sub

Figure 10-9. Fourth-order Runge-Kutta custom function for systems of differential equations, (folder 'Chapter 10 Examples', workbook 'ODE Examples', module 'RungeKutta3')

Figures 10-10, 10-11 and 10-12 illustrate the use of Runge3 to simulate some complex chemical reaction schemes. Figure 10-10 shows concentration vs. time for the consecutive first-order reaction scheme

A —> B—» C for which the differential equations are time, seconds

Figure 10-10. Runge-Kutta simulation of consecutive firstlorder reactions, (folder 'Chapter 10 Examples', workbook 'ODE Examples', worksheet 'A->B->C')

time, seconds

Figure 10-10. Runge-Kutta simulation of consecutive firstlorder reactions, (folder 'Chapter 10 Examples', workbook 'ODE Examples', worksheet 'A->B->C')

The parameters used in the simulation were [A]0 = 5.00 x 10~3 mol L~', k = 0.5 s-1 and k2 = 0.4 s"1.

Part of the spreadsheet is shown in Figure 10-11. The formulas for the derivatives, in cells G10, H10 and 110, are

and the formulas in cells J11, K11 and L11 are =Runge3(A10,J10:L10,G10:l10,A11-A10.1) =Runge3(A10,J10:L10,G10:l10,A11-A10,2)

 A 1 G i. H I i I J K L 8 9 time Using Rimye-Kiitta custom function (sec) tl( AJ/ctf (i[B] (It d|C|dt [fl]t [BJt [CJt 10 0.0 -2.50E-03 2.50E-Q3 O.OOE+OO 5.00E-03 O.OOE+OO 0 00E+00 1.88E-05 11 0.2 -0.00226 0,00203 1 83E-04 4.52E-03 4.57E-04 12 0.6 -1.85E-03 1.39E-03 4.58E-04 3.70E-03 1.15E-03 1.51 E-04 13 1,0 -1.52E-03 8.78E-04 6.38E-04 3.03E-03 1.59E-03 3.73E-04 14 1.4 -1.24E-03 4.95E-04 7.46E-04 2.48E-03 1 87E-03 6.52E-04 15 1.8 -1 02E-03 2.15E-04 8.02E-04 2.03E-03 2.00E-03 9.63E-04 15 2,2 -8.32E-04 1.31E-05 8.19E-04 1.66E-03 2.Û5E-03 1.29E-03 17 2.6 -6.81 E-04 -1.28E-04 8.09E-04 1.36E-03 2.02E-03 1.61 E-03 1B 3.0 -5.58E-04 -2.23E-04 7.81 E-04 1.12E-03 1.95E-03 1.93E-03 19 4.0 -3.38E-04 -3.27E-04 6.S5E-04 6.77E-04 1 66E-03 2.66E-03 20 5.0 -2.05E-04 -3.27E-04 5.32E-04 4.11 E-04 1 33E-03 3.26E-03 21 6.0 -1.25E-04 -2.84E-04 4.09E-04 2.49E-04 1 02E-03 3.73E-03 22 7.0 -7.56E-05 -2.30E-04 3.06E-04 1.51 E-04 7.65E-04 4.08E-03 23 8.0 -4.59E-05 -1.78E-04 2.24E-04 9.18E-05 5.61 E-04 4.35E-03 24 9.0 -2.78E-05 -1.34E-04 1.82E-04 5.57E-05 4 05E-04 4.54E-03 25 10.0 -1.69E-05 -9.89E-05 1,16E-04 3.38E-05 2.89E-04 4.68E-03

Figure 10-11. Spreadsheet for the Runge-Kutta simulation of consecutive first-order reactions, (folder 'Chapter 10 Examples', workbook 'ODE Examples', worksheet 'A->B->C')

Figure 10-11. Spreadsheet for the Runge-Kutta simulation of consecutive first-order reactions, (folder 'Chapter 10 Examples', workbook 'ODE Examples', worksheet 'A->B->C')

The arguments of the function can be entered in other ways. Two of these are illustrated in rows 12 and 13 of the spreadsheet. If the derivatives are located in non-adjacent cells, the deriv_formulas argument can be entered as a non-adjacent selection, as illustrated by the formula in cell J12:

The cell references must be enclosed in parentheses and separated by commas. The function can also be entered as an array formula, as in cells J13:L13

{=Runge3(A12,J12:L12,G12:l12,A13-A12)}

In this simulation, the largest errors are about 0.05%.

Figure 10-12 shows a second example, concentration vs. time for a second-order autocatalytic reaction scheme. An autocatalytic reaction is one in which a product acts as a catalyst for the reaction. The reaction has two pathways: an uncatalyzed path (A-»B) and an autocatalytic path (A + B —> 2B). The rate law (the differential equation) is

A ]Jdt = d[B]Jdt = &o[A], + kx [A],[B], ( 10-27)

The parameters used in the calculation were: ko = 1.00 x 10"4 s"1, kx = 0.50 M"' s"1, C = 0.0200 M. The spreadsheet can be examined on the CD-ROM.

0.020

0.010

0.000

0 200 400 600 800 1000 1200 time, seconds

Figure 10-12. Runge-Kutta simulation of second-order autocatalytic reaction, (folder 'Chapter 10 Examples', workbook 'ODE Examples', worksheet 'Autocatalytic') 0 200 400 600 800 1000 1200 time, seconds

Figure 10-12. Runge-Kutta simulation of second-order autocatalytic reaction, (folder 'Chapter 10 Examples', workbook 'ODE Examples', worksheet 'Autocatalytic')