Problem
You've seen how to apply both the RungeKutta method and Euler's basic method using VBA for first and secondorder differential equations, but now you need to solve a system of coupled equations . You're wondering how to extend what you've learned so far to deal with this more complicated problem.
Solution
You can apply the same VBA techniques discussed earlier with some modifications to deal with multiple equations.
Discussion
Consider this system of coupled nonlinear differential equations with initial conditions:
Pretty, huh? This problem originally consisted of two coupled secondorder equations that were reduced to four firstorder equations using the same technique discussed in Recipe 11.2. This is a common technique for reducing the order of differential equations, making them more amenable to solving. This particular problem was adapted from a problem that appears in An Introduction to Linear and Nonlinear Finite Element Analysis, by Prem K. Kythe and Dongming Wei (Birkhauser). The problem represents a finite element solution of the nonlinear vibrations of a metal rod. Since this is a nonlinear problem, it requires a numerical solution. In their book, Kythe and Wei present some Matlab code to solve this problem using Euler's method .
Euler's method
I implemented Euler's method in VBA to illustrate how to solve this same system in Excel. Example 114 shows my VBA code.
Example 114. DoEuler for nonlinear system
Public Sub DoEuler( ) Dim y(1 To 4) As Double Dim t(1 To 4) As Double Dim n As Integer Dim dt As Double Dim time As Double y(1) = 29.8613 y(2) = 0 y(3) = 30.0616 y(4) = 0 n = 5000 dt = 5 / n time = 0 For i = 1 To n t(1) = y(1) + dt * dy1dt(y(1), y(2), y(3), y(4)) t(2) = y(2) + dt * dy2dt(y(1), y(2), y(3), y(4)) t(3) = y(3) + dt * dy3dt(y(1), y(2), y(3), y(4)) t(4) = y(4) + dt * dy4dt(y(1), y(2), y(3), y(4)) y(1) = t(1) y(2) = t(2) y(3) = t(3) y(4) = t(4) time = time + dt ActiveSheet.Cells(i + 1, 1) = time ActiveSheet.Cells(i + 1, 2) = y(1) ActiveSheet.Cells(i + 1, 3) = y(2) ActiveSheet.Cells(i + 1, 4) = y(3) ActiveSheet.Cells(i + 1, 5) = y(4) Next i End Sub 
As you can see, this is not a particularly complicated subroutine. It's very similar to the code in Example 112 that I showed you for solving a firstorder equation using Euler's method. The major difference in this case is that we have four coupled equations instead of a single equation.

I have to confess that there's a little more to this code than is shown in Example 114. Notice the function callsdy1dt, dy2dt, dy3dt, dy4dtin the lines of code that implement Euler's integration formula (the lines just after the For statement). These functions compute the derivatives for each y. The derivative expressions are taken directly from the system of equations shown earlier. The corresponding VBA functions are shown in Example 115.
Example 115. Derivative functions
Public Function dy1dt(y1 As Double, y2 As Double, y3 As Double, y4 As Double) As Double dy1dt = y2 End Function Public Function dy2dt(y1 As Double, y2 As Double, y3 As Double, y4 As Double) As Double dy2dt = (15600# * (2# * ((Abs(y1)) ^ 0.74) + 3# * ((Abs(y3  y1))^0.74))) _ / 7# * y1 + (3# * 15600# * ((Abs(y3  y1)) ^ 0.74)) / 7# * y3 End Function Public Function dy3dt(y1 As Double, y2 As Double, y3 As Double, y4 As Double) As Double dy3dt = y4 End Function Public Function dy4dt(y1 As Double, y2 As Double, y3 As Double, y4 As Double) As Double dy4dt = (15600# * (((Abs(y1)) ^ 0.74)  5# * ((Abs(y3  y1)) ^ 0.74))) _ / 7# * y1  (5# * 15600# * ((Abs(y3  y1)) ^ 0.74)) / 7# * y3 End Function 
This code runs just fine. The problem is that this is a particularly troublesome system of equations for Euler's method. The solution diverges dramatically as it marches along in time. Check out Figure 115 to see what I mean.
The chart shows y1 and y3the displacements of the rod at its center and free endover time. The results should look sinusoidal. However, you can clearly see that this solution is not a sinusoid. This is wonderful example of the sort of problem where Euler's method, though easy to implement, is a poor choice.
RungeKutta method
The RungeKutta method presents a better alternative to Euler's method for this problem. In Recipe 11.2, I showed you how to reduce a secondorder equation to two firstorder equations and then use the RungeKutta method to solve one and use a Euler step to solve the other. In this new example, I'll show you how to apply the RungeKutta method for all four equations at each time step. To this end, I prepared a new subroutine, DoRungeKutta, shown in Example 116. DoRungeKutta also makes calls to the derivative functions shown earlier in Example 115.
Figure 115. Poor results using Euler's method
Example 116. RungeKutta for nonlinear system
Public Sub DoRungeKutta( ) ' Declar Local variables Dim y(1 To 4) As Double Dim k1(1 To 4) As Double Dim k2(1 To 4) As Double Dim k3(1 To 4) As Double Dim k4(1 To 4) As Double Dim n As Integer Dim dt As Double Dim time As Double Application.ScreenUpdating = False ' Initialize variables y(1) = 29.8613 y(2) = 0 y(3) = 30.0616 y(4) = 0 n = 5000 dt = 5 / n: time = 0 ' Perform iterations For i = 1 To n k1(1) = dt * dy1dt(y(1), y(2), y(3), y(4)) k1(2) = dt * dy2dt(y(1), y(2), y(3), y(4)) k1(3) = dt * dy3dt(y(1), y(2), y(3), y(4)) k1(4) = dt * dy4dt(y(1), y(2), y(3), y(4)) k2(1) = dt * dy1dt(y(1)+k1(1)/2#, y(2)+k1(2)/2#, y(3)+k1(3)/2#, _ y(4)+k1(4)/2#) k2(2) = dt * dy2dt(y(1)+k1(1)/2#, y(2)+k1(2)/2#, y(3)+k1(3)/2#, _ y(4)+k1(4)/2#) k2(3) = dt * dy3dt(y(1)+k1(1)/2#, y(2)+k1(2)/2#, y(3)+k1(3)/2#, _ y(4)+k1(4)/2#) k2(4) = dt * dy4dt(y(1)+k1(1)/2#, y(2)+k1(2)/2#, y(3)+k1(3)/2#, _ y(4)+k1(4)/2#) k3(1) = dt * dy1dt(y(1)+k2(1)/2#, y(2)+k2(2)/2#, y(3)+k2(3)/2#, _ y(4)+k2(4)/2#) k3(2) = dt * dy2dt(y(1)+k2(1)/2#, y(2)+k2(2)/2#, y(3)+k2(3)/2#, _ y(4)+k2(4)/2#) k3(3) = dt * dy3dt(y(1)+k2(1)/2#, y(2)+k2(2)/2#, y(3)+k2(3)/2#, _ y(4)+k2(4)/2#) k3(4) = dt * dy4dt(y(1)+k2(1)/2#, y(2)+k2(2)/2#, y(3)+k2(3)/2#, _ y(4)+k2(4)/2#) k4(1) = dt * dy1dt(y(1) + k3(1), y(2) + k3(2), y(3) + k3(3), y(4) + k3(4)) k4(2) = dt * dy2dt(y(1) + k3(1), y(2) + k3(2), y(3) + k3(3), y(4) + k3(4)) k4(3) = dt * dy3dt(y(1) + k3(1), y(2) + k3(2), y(3) + k3(3), y(4) + k3(4)) k4(4) = dt * dy4dt(y(1) + k3(1), y(2) + k3(2), y(3) + k3(3), y(4) + k3(4)) y(1) = y(1) + (k1(1) + 2 * k2(1) + 2 * k3(1) + k4(1)) / 6 y(2) = y(2) + (k1(2) + 2 * k2(2) + 2 * k3(2) + k4(2)) / 6 y(3) = y(3) + (k1(3) + 2 * k2(3) + 2 * k3(3) + k4(3)) / 6 y(4) = y(4) + (k1(4) + 2 * k2(4) + 2 * k3(4) + k4(4)) / 6 time = time + dt ActiveSheet.Cells(i + 1, 1) = time ActiveSheet.Cells(i + 1, 2) = y(1) ActiveSheet.Cells(i + 1, 3) = y(2) ActiveSheet.Cells(i + 1, 4) = y(3) ActiveSheet.Cells(i + 1, 5) = y(4) Next i Application.ScreenUpdating = True Application.ActiveWorkbook.Save End Sub 
The line Dim y(1 To 4) As Double declares a fourelement array variable in which to store the computed ys at every time step. The next several declarations consist of array declarations for the RungeKutta kvalues. Each of the four equations requires four empty kvalues.
After setting the initial y conditions and initializing other utility variables, the subroutine enters a For loop. The guts of the For loop contain a straightforward application of the RungeKutta formulas for each equation in the system being solved. As in earlier examples, output is returned to the active spreadsheet at each time step.
Running DoRungeKutta gives much better results than the earlier attempt that used Euler's method. Figure 116 shows the new results.
Figure 116. Good results using the RungeKutta method
These results look far better than those obtained earlier. Clearly, the results shown in Figure 116 for y1 and y3 are sinusoidal as expected. The difference between these results and those obtained using Euler's method are dramatic and effectively illustrate the relative accuracy of higherorder methods, such as the RungeKutta method, over Euler's basic method.

Using Excel
Getting Acquainted with Visual Basic for Applications
Collecting and Cleaning Up Data
Charting
Statistical Analysis
Time Series Analysis
Mathematical Functions
Curve Fitting and Regression
Solving Equations
Numerical Integration and Differentiation
Solving Ordinary Differential Equations
Solving Partial Differential Equations
Performing Optimization Analyses in Excel
Introduction to Financial Calculations
Index