Numerical Differential Equations


Introduction

Often when doing transient analysis of a system, being able to quickly generate a simulation of first or second order differential equations (that are not necessary linear) is very useful. This tutorial comes with an example spreadsheet and the steps by which you can easily make a simulation within Microsoft Excel.

Model Formulation

Suppose that the differential equation is non-linear and second order, such as in the case of a body moving through fluid. In this example, we will consider a vehicle using a DC motor driving in a straight line with fluid drag. The goal is to find out how long it will take for the car to get to maximum velocity, and how much distance this process would need.

First of all a simple force analysis gives us the equation below

$$ F_{motor}-F_{drag} = mx'' $$

For steady state speeds, DC motors exhibit a linear decline in torque with angular velocity at an applied voltage. This translates to linear force and velocity by a scaling factor. Thus, a simple model of the motor force, which neglects the transient behaviour of the motor itself, is given below.

$$F_{motor} = a_0-a_1x' $$

Drag is simplified to a constant drag coefficient multiplied by \(x'^2\).

$$F_{drag} = C_dx'^2 $$

Therefore, by substitution

$$ a_0-a_1x'-C_dx'^2 = mx'' $$

The equation is formulated and will be solved in the next section via Euler's method.

Implementation

Euler's method is a first order method to approximate differential equations and therefore needs a relatively small time step in order to be accurate. There are certainly better ways to numerically solve differential equations however, it is extremely easy to apply Euler's method in Excel for one dimensional problems.

First as with any differential equation, either boundary conditions or initial values are necessary to solve the equation numerically. In this case, let us say that \(x = 0\) and \(x'=0\) define the initial conditions.

The equations below come from integration of a constant acceleration to find the velocity and displacement. This is why the time step must be small, otherwise the assumption will produce inaccurate results.

$$ x_0 = 0 $$

$$ x'_0 = 0 $$

$$ x''_n = \frac{1}{m} \left(a_0-a_1x'-C_dx'^2\right) $$

$$ x'_n = x''_{n-1}\cdot dt + x'_{n-1} $$

$$ x_n = x''_{n-1}\cdot dt^2 + x'_{n-1}\cdot dt + x_{n-1} $$

$$ t_n = t_{n-1} + dt $$

Now all we have to do is place \(t,x,x',\) and \(x''\) into separate columns in a table, specify the initial conditions, and input the formulas. An example of such a table is shown below. I also included the table of constants on the left. Note that you only need to type the formulas in the first and second rows, then you copy paste the second row downwards to extend the simulation.

Graph your results and that's it! With my parameters used above the graph for the velocity and displacement is shown on the right.

From the graph we can see that the velocity \(x'\) reaches an approximate maximum at 0.7 after 2.2 time increments have passed. The displacement needed for the velocity to reach a maximum is found to be approximately 1.2.

You can download this specific simulation spreadsheet below and alter the values to see the different results.