Ordinary Differential Equations (ODE) tool

Building a model

ModelRisk’s Ordinary Differential Equation (ODE) tool will numerically evaluate one or more variables over time that follow one or more ordinary differential equations . It allows the user to define any set of first order differential equations that can be described with Excel functions. One or more time stamps (specific points in time) can be specified for the evaluation of the variable(s). The interface will plot any variable against time or any two variables together. The user can also specify shock to the system at specific points in time that change a variable or model parameter.

Example

Click here to open the model for this example. A missile is fired from a height h and an angle θ and initial speed . The missiles speed reduces due to friction by an amount * where Friction is a friction constant and is its speed at time t. Model the speed and location of the missile during flight.

The diagram below shows this system. If you studied maths at college, you’ll probably remember that the easiest way to solve this problem is to consider the movement in the horizontal (x) and vertical (y) directions separately, since gravity only works in the vertical direction.

We can easily write down the differential equations for the horizontal and vertical speeds and locations (x,y) as follows:

where g is gravitational acceleration.

To use the ODE tool we need to translate each of these equations into text using standard Excel functions and notation:

We will also use the following initial values (using SI units - meters, seconds):

The ODE tool can be accessed through the ‘More Tools’ list in the ModelRisk ribbon:

This opens the ODE Tool interface:

The interface has three sections:

Left side: Entry of time horizons, equations, initial values for variables and (more advanced) any perturbations to the variable values during the time horizon.

Center: Graphs of the variables being modelled.

Right side: Statistics of the variables being displayed in the graph.

We will begin by entering the Start time (=0), Step (the size of the discrete time step used ∆t to approximate continuous time, = 0.01 seconds), and the ‘Time stamps’ (the times at which we would like ModelRisk to report the variable values. In this case we will just use one time stamp = 6 seconds:

Now we enter the differential equations by one of two means:

Clicking the ‘Add’ button and then directly typing each equation into the interface; or
Placing the equations in Excel and then clicking the ‘Add’ button and then the ‘Ref’ button four times to select the differential equations in cells B3:B6:

Next we enter the initial values of parameters, by typing the name of the variable in the left column and either again using ‘Add’ and ‘Ref’ buttons for the values, or typing those values directly into the interface:

If values are linked to spreadsheet cells, the variables can also of course be random samples from distributions. In this simple example, there are no discrete perturbations to the variables, so the last table is left blank.

All the data have been entered and the ‘Errors’ message box is empty, showing that all variables have been defined and the differential equations are syntactically correct. We can now save the model by clicking the ‘OK’ button and following the instructions on where ModelRisk is to place the different model components.

ModelRisk will take all the inputs and arrange them together into your spreadsheet, together with any links you created. It has also inserted a VoseODE array function into the model, shown highlighted here:

Clicking the ‘View function’ icon while one of the cells covered by the VoseODE function is active will bring back up the ODE dialog box:

The two variables and have been selected from the differential dialog using tick boxes and are therefore displayed graphically. One or two variables can be displayed at the same time and against a different x-axis variable selected at the bottom of the screen. In the following plot, the vertical distance y is plotted against the horizontal distance x. This shows the trajectory of the missile:

The statistics show that the missile reaches a maximum height of 30.98 meters above ground, and lands about 46 meters from where it was fired (where it crosses the x-axis).

Saving the model

You can save the model you have created by simply saving the Excel file that you have built. However, ModelRisk also allows you to save the set of differential equations in its own library. This will allow you to pull the set of equations up again when you have anew model to build. To save the equation set to the ODE library, enter a name for the model in the Name dialog box:

Click the Description button to add some descriptive text to explain the model:

Click OK, and then the Save button. When you next open the ODE tool interface you can click the Load button and select the template from the library:

The ODE tool will then load the entire model you built, including time stamps and initial values, but without links to a spreadsheet:

Advanced feature: perturbations

ModelRisk allows the user to specify that a change can occur to one or more variables at a particular point in time. For example, the set of ODE equations could be a PK/PD (Pharmacokinetic / Pharmacodynamic) model describing the effect of a particular regime of drug application. The perturbations would then represent sudden changes in the concentration of the drug due to taking a pill, or receiving an injection.

For simplicity, we will stick with the missile model. Imagine that after 1 second of flight the non-slip coating on the missile has worn off and the Friction parameter now changes to a value of 0.95 /s. We enter this in the ODE interface as follows:

The effect on the missile path is obviously going to be rather dramatic. In terms of horizontal speed it will slow down very fast, and in vertical speed gravity will dominate. This can be seen by looking at different graphical plots in the ODE interface:

The trajectory of the missile will show a much steeper drop-off (compare with other plots above):

Using the ODE Tool

Using the ODE tool to find a target time

We can find out when and where the missile hits the ground by using Excel’s Solver:

This is determining the time stamp that gives y=0:

The model calculates that the missile will land 46.29 metres away after 4.62 seconds.

Using the ODE Tool to fit to data

Click here to open the model for this example. Imagine that we did not know the Friction constant, but were able to track the position of the missile producing the following data and wished to estimate the value of Friction from these data:

We can set up the model with four time stamps (t={1,2,3,4}):

and click ‘OK’ to enter the function into the spreadsheet:

In this version of the model cell J24 calculates the sum of the squared errors between the observed and estimated x and y values. A dummy Friction value of 0.3 is entered into the model, and we can use Excel’s Solver to vary Friction until the smallest Total error is achieved:

The Solver returns the value of 0.100015, which is almost precisely the 0.1 value for Friction that was used to generate the data used in this example.

Using the ODE tool with VBA macros and Solver

Click here to open the model for this example. Imagine that we wish to estimate the variation in range that the missile will travel as a result of small variations in initial muzzle velocity and the coefficient Friction. We will replace the values in the spreadsheet with the following distributions:

: =VoseLognormal(25,0.7)

Friction: =VoseLognormal(0.1,0.029)

Since the model requires Solver to run, recalculating the spreadsheet several times to arrive at a solution, this model has a VBA macro that performs the following tasks:

Recalculate the spreadsheet;
Copy random sample values from the above distributions into the cells for and Friction;
Using these copied values, run Solver to find the value of FlightTime that gives a y value of zero

The macro is run before each simulation sample of the model by specifying it in the ModelRisk Simulation Settings dialog box:

The Range (x) is stored as an output and the sampled values for and Friction are stored as inputs. After 5000 samples, we get the following output result:

Simulation results for the Range show it will lie between about 35m and 60m. A sensitivity plot verifies that the larger the friction coefficient, the smaller the range (negative correlation) and the higher the muzzle velocity the greater the range (positive correlation) which gives us some intuitive confidence that the model is behaving correctly.