This article focuses on two types of integrators for simulating the motion of an object in Matlab (or Octave). Videos and downloadable source code are at the end of this article. Note that a free alternative to Matlab, called Octave, can also be used to run the software and this is covered here and in one of the videos.
If you’re taking on the task of building a simulation from scratch, the first place you start is in deciding what integrator you’ll be using for propagating the trajectory – whether it’s for a simple case (such as this one) or for a complex aerospace system traversing the atmosphere and space itself. The next step – required by the integrator – is modeling the forces that are acting upon the object that you’re simulating. Once you’ve accomplished these objectives, you’ve created the core engine that drives your simulation.
To keep it simple, we’ll look at a 1-Degree-of-Freedom (DOF) system – in this case a spring / damper system, that moves along the X-axis, and is acted upon by an external, time-dependent forcing function. The spring produces a force based on the displacement from its natural state (positive or negative, depending on whether it’s compressed or extended). The damper produces a force based on the velocity of the object (the sign, negative or positive, is determined by the direction of the velocity).
While we’ll be working with Matlab (and demonstrating with Octave as well), one of the integrators can be used in any software development language.
Matlab has several available integrator functions – the most popular is ODE45 (ODE = Ordinary Differential Equations) which completes the entire trajectory simulation process in one call to the function. An example of its use (from the actual code base in this project) is shown below. It’s very elegant but simple in that everything is accomplished in one call to this function – the entire trajectory set (position and velocity, as well as time profile) is returned in the stateVec array and time vector on the left hand side.
Runge-Kutta 4th Order
The other integrator is the Runge-Kutta 4th order integrator (also very common and popular) – in this case the trajectory is propagated one time-step at a time. As shown in the code example below, the integrator is called for each time step in the simulation and thus the solution is propagated forward with each call to the function. This is not a Matlab or Octave native function and thus has to be hand-coded (or the basic algorithm copied from another source and converted into Matlab’s scripting format).
Pros and Cons for Each Approach
The advantage with this approach is that it’s a single call to the function (no iterative calling of the function throughout the entire trajectory as is the case with the Runge-Kutta approach), and it has a variable time-step capability such that it will determine the best time step for each phase of the trajectory. This latter feature frees up the user from trying to figure out what time step size should be used to generate accurate simulation results.
The disadvantage is portability – that is that if you decide to convert the code base to say Java or C++, the Matlab ODE45 function cannot be carried over. Thus you’ll have to either use a built-in Java or C++ integrator function or else use a hand-coded method such as Runge-Kutta.
Runge-Kutta 4th Order
The advantage with this popular and time-tested approach is portability – it can easily be converted in any language. The disadvantage is that it is a fixed time step method which means that you must have a reasonable idea of the value of the time step.
The two key elements are: 1) the integrator, and 2) deriving the force equations acting on the model. The illustration below shows the system diagram and the associated function calls for the spring force (kForce), the damper force (dForce), and the forcing function (fForce). The integrator requires the acceleration values so we simply divide the forces by the mass of the object to obtain the acceleration values.
The forces and acceleration are simply modeled as shown in the above diagram. The previous state (position and velocity) is brought in via the stateVec array and loaded to x and xd. The spring force, kForce, is computed using the spring constant, k, and the displacement, x. The damper force is computed using the damper constant, c, and the velocity, xd.
Note that we have to have reference signs for the accelerations. If the object moves along the positive X-axis, then the spring is compressed and produces a force in the negative direction of the X-axis. The damper force is dependent on the velocity (magnitude drives the force value, the direction determines the sign).
Results in Matlab
Identical results can be obtained when using either integrator as shown below. In this case, the Matlab ODE45 integrator was run with a variable time step, while the Runge-Kutta integrator was run with a very small fixed time step of 0.01 seconds.
Position and Velocity Profiles
Given that ODE45 data is plotted in blue and the Runge-Kutta data is plotted in red, it’s hard to see any differences in the trajectories – they are basically identical.
Time Step Profiles
The time step profiles for each integrator are shown below. Note with the variable time step profile on the left (ODE45), the time step starts out very small because the forces change rapidly in the beginning (due to the relatively high frequency of the forcing function). As the forces change less over time, the time step size is increased.
Of course, the Runge-Kutta time step (profile shown on the right) is fixed so it never changes regardless of the state of the forces acting on the object.
Results in Octave
As with Matlab, identical results can be obtained when using either integrator as shown below. In this case, the Octave ODE45 integrator was run with a variable time step, while the Runge-Kutta integrator was run with a very small fixed time step.
Position and Velocity Profiles
The results are similar to those from the Matlab test runs – it’s hard to see any significant differences in the trajectories – they are basically identical.
Octave ODE45 Function Call
One change had to be made to the ODE45 integrator implementation – that is that the tolerances had to be lowered because it appears that the default values are too high for the example test case (the high frequency of the forcing function requires tighter tolerances in the variable time step algorithm). Thus an ODE45Wrapper function was made to contain the ODE45 function calls for both Matlab and Octave. The difference is that the tolerances are specified prior to the ODE45 call, in Octave, as shown below.
Time Step Profiles
The time step profile is different than that of Matlab but the trend is the same – it starts small and gets larger over time as the frequency of the force profiles decrease.
For more detailed information, you can watch any or all of the videos below. And as previously mentioned, the complete code base is available for download after the videos sections.
Video 1 – Running the Integrators in Matlab
This video below shows how to run the code when using Matlab. Click on the lower right square icon (next to the sound / speaker icon) to enlarge the video to almost the size of the monitor in order to more easily view it.
Video 2 – Running the Integrators in Octave
This video below shows how to run the code when using Octave. Click on the lower right square icon (next to the sound / speaker icon) to enlarge the video to almost the size of the monitor in order to more easily view it.
Video 3 – Code Walk-Through
The video below is a code “walk-through” if you’re interested in getting some details of the various functions. Click on the lower right square icon (next to the sound / speaker icon) to enlarge the video to almost the size of the monitor in order to more easily view it.
The code base for this article (directories and source code) is available for download as a zip file from the link below. Feel free to use this code for your own purposes with no obligation whatsoever to me. However, if you feel that it’s been beneficial to your efforts, please refer friends and colleagues to this blog when appropriate – thank you!!
If you have questions or comments, please leave comments – or you can email me directly at email@example.com