Solving ODEs on Linux
Many problems in science and engineering are modeled through ordinary differential equations (ODEs). An ODE is an equation that contains a function of one independent variable and its derivatives. This means that practically any system that changes over time can be modeled with an ODE, from celestial mechanics to chemistry reaction rates to ecology and population modeling.
Because of this ubiquity, many tools have been developed through the years to help solve and analyze ODEs. In this article, I take a look at one of the tools available on Linux: Model Builder. The project is hosted on SourceForge, so you always can build it from source, but most distributions should have a package available. On Debian-based distros, you can install it with the command:
sudo apt-get install model-builder
It also installs several Python modules to support the tasks it can handle. If you do decide to build from source, you will need to handle these dependencies yourself.
Included with the source is a directory of examples. You can use them as a starting point and to gain some ideas of what you can do with Model Builder. Documentation is a bit sparse, so you may need to get your hands a little dirty to take the most advantage of what is possible with Model Builder.
To start Model Builder, you either can click on its menu item in your
desktop environment or run the command PyMB
from a terminal
window. When the main window pops up, you are presented with a template
where you can define the problem you are analyzing (Figure 1). The main
pane, titled Differential Equations, is where you can define the set
of ordinary differential equations that you are trying to solve. The
general form of these equations is dy/dt = f(y,t).
Figure 1. When Model Builder starts, you can set several parameters and the equations you want to analyze.
If your system depends on different levels of differentiating the dependent variable, you always can rewrite it as a system of ODEs. When you give Model Builder your system, you need to write out only the right-hand side of the above equation. This equation can contain essentially any function or expression that NumPy understands, since Model Builder uses Python to do the heavy lifting.
Because Model Builder is designed to handle systems of equations, you need to define the y portion as elements of a list. So the y variable for the first equation is labeled as y[0]; the y variable for the second equation is labeled y[1] and so on. These are called the state variables.
The pane to the right of the equation window is where you can place any parameters that you need, one per line. They can be used in the equation window, where they are labeled as p[0], p[1] and so on. If you want to use time in either the parameters or equations that you have defined, you just need to use the t variable.
Because Python is used in the back end, you even can use lambda functions to define more complex structures. You may want to take a look at the documentation available on the NumPy site to see what options are available.
Below these two panes is where you define the rest of the options for your problem. In the Initial values box, you can enter the initial values for each state variable at the time t=0. They need to be separated with a space and put in the order of the equations given in the equation pane.
Below the Initial values, you can enter the start time, the end time and the time step to use in the solution. The critical time steps box is usually left empty, so let's leave it alone here. The first step box is the size of the first step. Usually, you should leave this as 0 to allow for automatic determination. The minimum and maximum step size boxes set these variables that are used in the variable step size algorithm. Typically, you should leave these as 0 as well to allow for automatic determination. The full output check box will print out more useful information about the integration in the results spreadsheet.
Once everything is entered, all you need to do is click the Start icon, and the integration will be calculated. If this is a system that you will want to work with over time, you can click on the menu item File→Save to save the model to a file. This file format is an XML file, so you could edit it with a text editor if you want. When you are ready to do more work with it, you can load it by clicking on File→Open.
Once the calculations are done, which may be fast for simple problems, a results window will pop up (Figure 2). matplotlib handles this graph window, so you can manipulate it just like any other matplotlib window. This includes panning, zooming or changing the plot window. You also can save the resulting plot as an image file in one of several different formats.
Figure 2. Once you finish defining the problem and run the integration, a result window pops up with a graph of the integration.
Going back to the main window, let's look at some other available tools. Clicking on the Show equations icon pops up a window where you can see the equations typeset (Figure 3). Beside this icon is the Results icon. Clicking on that pops up a spreadsheet of all of the results from your integration (Figure 4). The columns of data include the time, the value of y[0] and the step sizes, among other things. You can select a couple columns by holding down the Ctrl key and clicking on the column headers. Then, click on the plot button to plot them in a new window. You can get a power spectrum for any one column by selecting one of interest and clicking on the Spectrum icon. This pops up two new windows, the first a power spectrum of the column (Figure 5) and the second a spectrogram of the column (Figure 6).
Figure 3. You always can get a typeset display of your equations to verify what they should look like.
Figure 4. You can pull up all of the results of your integration and do further analysis.
Figure 5. You can generate a power spectrum of any column of your results.
Figure 6. You also can generate a spectrogram of your results.
The last tool available is a wavelet transform. When you select a column, you can apply a continuous wavelet transform to the data. When you are done with Model Builder, you can save this data into a comma-separated values (CSV) file from the spreadsheet window. Then, you can import it into other tools, like R, to do even further analysis.
Now that you have seen the options available in Model Builder, hopefully you will consider it when looking at ODE problems. It provides a pretty simple interface to the tools available in Python to solve ODEs numerically. Although other more powerful tools are available, Model Builder fits into the niche of experimenting quickly with different equations and playing with ideas.