Problem
You need to numerically solve a secondorder differential equation of the form:
Solution
This is a standard initial value problem, and you can implement any of a number of standard numerical integration techniques to solve it using Excel and VBA.
Discussion
Consider this equation and initial conditions:
Physically, this equation represents the equation of motion for an object moving under some applied thrust, T.[*] Here, m represents the mass of the object, Cd represents a drag factor, and s represents the displacement (position) of the object. This equation allows us to compute the motion of this object.
[*] This problem was adapted from an example I used in one of my previous books, Physics for Game Developers (O'Reilly).
To solve this equation, it helps to rewrite it in the form of two firstorder equations. Let v = ds/dt represent the velocity of the object. Now we have the following:
These equations represent two coupled firstorder equations, and you can apply any of a number of numerical integration schemes to solve them both.

I decided to use the RungeKutta method for this example. This method is based on taking more terms in the Taylor series expansion of a function, as I explained in Recipe 11.1. The caveat here is that you need to apply further Taylor series expansions to estimate the higherorder derivatives. The RungeKutta approach reduces the truncation error to something on the order of (dt)5, as opposed to (dt)2 in the case of Euler's basic method (dt is the step size). This means you can usually use larger step sizes while still maintaining accuracy. The tradeoff here is more computations per step. But ultimately, with the larger step size you should be able to reduce the overall computational burden using this method.
The RungeKutta integration formulas are a little more complicated than the simple formula we used with Euler's basic method. The general equations are as follows:
y' represents dy/dx in these equations. You can see that every step requires four intermediate computations prior to predicting the next yvalue, that is, y(x + dx).
Example 113 shows the VBA subroutine I wrote implementing these formulas.
Example 113. Implementation of RungeKutta method
Public Sub DoRK2ndOrder( ) Dim T As Double ' Thrust Dim Cd As Double ' Drag coefficient Dim M As Double ' Mass Dim dt As Double ' Time step size Dim F As Double ' Force Dim A As Double ' Acceleration Dim Vn As Double ' Velocity at time t Dim Vn1 As Double ' Velocity at time t + dt Dim Sn As Double ' Displacement at time t Dim Sn1 As Double ' Displacement at time t + dt Dim time As Double ' Total time Dim k1 As Double ' RungeKutta k1 Dim k2 As Double ' RungeKutta k2 Dim k3 As Double ' RungeKutta k3 Dim k4 As Double ' RungeKutta k4 Dim n As Integer ' Counter controlling total number of time steps Dim C As Integer ' Counter controlling output of results to spreadsheet Dim k As Integer ' Counter controlling output row Dim r As Integer ' Number of output rows ' Extract given data from the active spreadsheet: With ActiveSheet dt = .Range("dt") T = .Range("T") M = .Range("M") Cd = .Range("Cd") n = .Range("n") r = .Range("r_") End With ' Initialize variables: k = 1 time = 0 C = n / r Vn = 0 Sn = 0 ' Start iterations: For i = 1 To n ' Compute k1: F = (T  (Cd * Vn)) A = F / M k1 = dt * A ' Compute k2: F = (T  (Cd * (Vn + k1 / 2))) A = F / M k2 = dt * A ' Compute k3: F = (T  (Cd * (Vn + k2 / 2))) A = F / M k3 = dt * A ' Compute k4: F = (T  (Cd * (Vn + k3))) A = F / M k4 = dt * A ' Compute velocity at t + dt: Vn1 = Vn + (k1 + 2 * k2 + 2 * k3 + k4) / 6 ' Compute displacement at t + dt using Euler prediction: Sn1 = Sn + Vn1 * dt ' Update variables: time = time + dt Vn = Vn1 Sn = Sn1 ' Output results to the active spreadsheet: If C >= n / r Then ActiveSheet.Cells(k + 1, 1) = time ActiveSheet.Cells(k + 1, 2) = Sn ActiveSheet.Cells(k + 1, 3) = Vn k = k + 1 C = 0 Else C = C + 1 End If Next i End Sub 
This is a relatively long subroutine compared to the one shown in Example 112 implementing Euler's method, so I commented this one to help you navigate through it.
As usual, several local variables are declared at the beginning of the subroutine. There are a bunch of them here, so I labeled each one in Example 113 to indicate their respective purposes.
After the variables are declared, several of them are populated with data retrieved from the active spreadsheet. The spreadsheet is shown in Figure 114.
Figure 114. Solution using the RungeKutta method
Values retrieved from the spreadsheet include the timestep size, thrust, mass, drag coefficient, number of iterations, and number of output rows. The block of code in Example 113 containing the With statement performs the data retrieval. Notice that I use cell names to refer to the spreadsheet cells containing the data of interest.

The next little section of code simply initializes a few variables. This is where the initial values for displacement and velocity specified in the problem are set.
Next, the For loop iterates computing estimates for displacement and velocity at each step. Taking a closer look at the code within the For loop, you can see how I applied the RungeKutta formulas to compute vt+dt. To make things a little more intuitive (at least to me), I broke the computations of the kvalues into a few steps, using the analogy of Newton's second law of motion. F in these computations stands for the sum of forces. For this problem, the sum of forces is composed of the thrust force and the drag force. A then represents the object's acceleration, which is equal to the net force divided by the object's mass. Finally, the kvalues are computed by multiplying the acceleration by the time step.
I find it easier to think in terms of physical values like force and acceleration whenever possible. Knowing that acceleration is equal to dv/dt, you can readily see that acceleration corresponds to y' in the RungeKutta formulas for this particular problem. This is because we broke the original equation up and are dealing with the first one involving velocity.
After all of the kvalues are obtained, vt+dt is estimated using the RungeKutta formula discussed earlier. Next, a simple Euler estimate is made for the corresponding displacement, st+dt.
Finally, the results are sent to the active spreadsheet in the same manner I showed you in Recipe 11.1. Figure 114 shows the results, including a plot of velocity and displacement versus time, using a time step of 0.005. You should play around with this time step to see how the accuracy of the results is affected. (If this were a class, I'd assign that exercise for homework.)
I should also mention that I used a button control, as you can see in Figure 114, to call the subroutine when the button is clicked. Refer back to Recipe 9.3 to learn how to add this sort of button control and assign its functionality.
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