Recipe 11.3. Tackling Coupled EquationsProblemYou'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. SolutionYou can apply the same VBA techniques discussed earlier with some modifications to deal with multiple equations. DiscussionConsider 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 methodI 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
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
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 y_{1} and y_{3}the 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 methodThe 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 methodExample 116. RungeKutta for nonlinear system
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 methodThese results look far better than those obtained earlier. Clearly, the results shown in Figure 116 for y_{1} and y_{3} 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.

Recipe11.3.Tackling Coupled Equations
Similar products