Numerical integration is the process of utilizing numerical methods to estimate the value of a definite integral and more generally numerical solutions to differential equations. In many cases an equation may have an explicit solution that can be solved for, in which case you should just use that, but for many other equations there may not be and numerical integration can help solve them.
As mentioned before these integrators can see a lot of utility in scientific computation and simulation but I will be focusing on their usage for dynamic simulation today.
The math behind the three integrators I will be going over today is not that bad, in fact the three methods extend well on each other.
The Euler integrator is the simplest one that we will consider today. The basic principle is to repeat the following steps.
Which looks something like this:
Where is the amount to take a step with on each iteration (something you determine and often called dt). To be a bit more rigorous let's manually compute a few steps of the Euler integration for a function using .
This differential equation is one we can integrate explicitly so let's compare and see how we did.
Let's compare the values from the known solution to the results from Euler.
| Euler | ||
|---|---|---|
As we can see the results are not very impressive, we'll discuss more about the limitations of euler integration.
The Heun integrator extends on Euler by adding a second derivative estimate and taking the average of them. The general method is repeating the following steps:
This looks like:
Following the same function as prior let's manually compute two steps.
Let's go ahead and see how the Heun integrator stacks up against the Euler and true results.
| Euler | Heun | ||
|---|---|---|---|
As we can see it's still not doing great but it's closer than Euler.
The RK4 integrator extends on both Euler and Heun by computing 4 estimates of the required changes and taking a weighted average of them based on the butcher tableau.
This looks like:
Following the same function as prior let's manually compute two steps.
Let's go ahead and see how the RK4 integrator stacks up against the Euler, Heun, and true results.
| Euler | Heun | RK4 | ||
|---|---|---|---|---|
As we can see the RK4 did a much better job estimating the true value! The steps that we took were quite course for this particular equation but we'll see later on that the trend of the RK4 clearly beating the other two integrators continues into the implementation.
Let's apply three common numerical integrators to an undamped and unforced spring mass system. Let's first detail the differential equation for this system, it is as follows:
Where is the mass of the object, is the spring rate, is the term we desire to track the position of the mass, and is the acceleration of the mass. Let's start by defining this differential equation as a Python function that returns a matrix that contains the velocity and the acceleration .
With this structure we can do both integrations in one step! As we are starting with and we are trying to integrate up two power to just we need to integrate the acceleration to get velocity and integrate once more to get . With this matrix setup the matrix addition will update both terms at once!
Next we can implement the function that applies an Euler integration step to solve for the velocity and position of the mass. As a reminder this is the first of the integrators mentioned prior and only relies on the rate at that instant.
Which takes in a series of arguments and computes the estimated state one step ahead, notably the evaluate argument is where we can pass in the function, eval_derivs, for our equation derivative. This function only computes one step when called so we need to wrap it with an iterator to advance us through multiple steps. This can be done as follows:
This chunk of code sets up our initial conditions for our system ( and ), the step to advance over, the total time, and a matrix of size to store and over time. Although I mentioned that numerical integrators can be used to integrate functions with unknown explicit solutions this is not one of them so we can use the true solution to see how well our integrator works!
The result is not the best but we can lower the step to improve the accuracy.
Although reducing the time step helps the result is still pretty bad, but why is this the case? As it turns out the simplicity of Euler integration, while convenient, leaves a lot to be desired in terms of accuracy. The single estimate of the rates to predict further steps does not allow the integrator to handle concave/convex sections without inducing error. To better understand this let's take a look at 2 coarse integration steps using Euler on this function.
Consider the first arrow, at the value of and so integrating that with Euler we get the following:
So no matter what h is the very first step of this integration will never change the values of . This is a problem because we can see that the mass does start moving and no longer becomes zero right after . This is where the multiple terms of the Heun and RK4 integrator can help a lot!
We can follow the same steps we took for Euler and all we need to do now is implement heun in the code and pass that in as the integrator in our iterator, that looks like this:
And for the integrator all we need to change is what integrator function call is used.
The first term k1 is identical to Euler but we introduce k2 to get an estimate of the rates from evaluate at the result of the Euler step. We then take an average to estimate what the average rate from the start to end should be. This yields the following results:
This looks great compared to Euler! We still can see some error produced near the end but it's much better. Let's see what we can achieve with coarser and finer steps.
While the results are appreciably better than Euler we still see that we need to make the step quite small to get reasonable results. Going back to the lack of handling concave/convex regions lets look at the closeup.
The first step is nearly perfect but we can see that the second step overestimates the rates and so we over correct, let's see what the additional two terms of RK4 can do.
Similarly to Heun, all we need to do is implement rk4 and change the integrator in the iterator.
We start with k1 identical to Euler but then it gets interesting, as a reminder the terms can be interpreted as:
The rates are combined using a weighted average based on the butcher tableau. All we do now is update the integrator call in our iterator.
Let's see how well this performs!
Wow! The rk4 performs extremely well for the same time step as used prior, let's see what this looks like for other time steps.
Even at the coarsest time step we get a very usable result! This is great because for a coarser timestep we can save on compute (even if each step has more compute involved when compared to Euler/Heun).
This is no surprise given the prior plots but the close up shows the steps landing almost exactly on the true function.
With this example we will introduce a damping term and show how we would write each step but also dive a bit more into the error analysis of these methods. We will now be integrating the following:
First we need to redefine our eval_derivs function for this new differential equation.
Now we just need to call the integrators within the iterator as prior, we'll also be using the same initial condition. To refresh here is what this looks like again.
No changes need to be made to our integrators as they will operate the same as before.
Below are the results for Euler integration on this damped oscillator.
This is hardly a result and so this integrator is worthless at this timestep. The error analysis obviously agrees as visually it's apparent already.
Below are the results for Heun integration on this damped oscillator.
Heun does much better than Euler but seems to produce a bit of a phase shift or lag in the results. The error analysis shows this did much better but there is much more to be desired still (especially for percent error).
Below are the results for RK4 integration on this damped oscillator.
RK4 does wonderful here, we can not appreciate much error when looking at the points of the steps. The error analysis shows very small absolute and percent errors. Although this is the case I would still like to point out something interesting about the percent error.
It's easy to try to chase down the values for both percent and absolute error but I want to offer some advice on how to perceive errors and determine what is "real". Let's imagine we were tracking the position of a car in a grid using an integrator. Far from the origin (0,0) a small tracking error of 1 meter for example could be a very small error whereas the same tracking error near the origin could be huge!
Let's say the true x position of the car is 400m from the origin with a tracking error of 1m, let's compute the errors. Obviously the absolute error is 1m so let's just compute the percent error.
Let's consider the same tracking error but the x position is now 5m. The absolute error is the same in both cases but as we'll see the percent error is very different.
To close the error conversation, it helps to think about what type of error measurement is important. In something like a car where you're concerned about staying in a lane the absolute error is all that matters whereas the percent error could give you a grossly over or unexaggerated idea of how well you're doing. There are other cases where percent error can be more useful such as when the values do not approach a sensitivity point near 0 or when relative comparisons are to be made.
This is merely an introduction to using these methods but if there is interest we can extend this material to cover more points such as:
Numerical integrators can be fast and flexible solutions to solving difficult ODEs in a programming environment. Although they seem simple in principle much of the difficulty can be in validating the performance of the integrators as we've seen simple integrators struggle with functions with continuously changing derivatives.
I hope you enjoyed and learned something new! Let me know if you have any questions or suggestions for improvements.