Continuous-time Dynamical Systems Runge-Kutta Numerical Integration
Dr Computer Science University of Auckland
Euler error
Copyright By PowCoder代写 加微信 powcoder
We have covered the Euler method for numerically integrating differential equations.
We talked a bit about how this method is not very accurate, but that we can increase its accuracy by decreasing the ‘step-size’ (¦¤t).
If we take a function for which we know the exact integral (via algebraic methods), we can visualize this increase in accuracy.
Euler error
The smaller the step-size, the more calculations you have to perform to simulate a system for a given duration.
# iterations = duration / stepsize
Which implies a practical lower-limit to ¦¤t (due to computational cost).
It’s also worth pointing out that error in these systems is generally accumulated, so the longer the duration, the less accurate the (later parts of the) simulation.
Other methods
A variety of other methods exist for numerically integration. One of the most frequently used (and easier to understand) is known as the Runge-Kutta method.
Actually, RK refers to a whole class of methods, but there is one specific version of it, known as RK4, or just the “classic Runge-Kutta” method that we’ll focus on.
As we shall see, it uses more calculations per time step, but allows for larger time steps, ultimately reducing the total computational cost necessary for simulating a reasonably accurate trajectory.
RK4 works by estimating how the derivative f(y,t) changes on its way to the next moment, i.e. as yt goes to yt+¦¤t.
Let’s take a look at how it works. To do so, we walk through one iteration of RK4, where we are transitioning from t = 0, to t = ¦¤t. We will use a large ¦¤t=1, as picking a large ¦¤t exaggerates the visualization making things easier to see.
1. evaluate k1, the derivative at t0. 2. evaluate k2, the derivative at t¦¤t/2
(as estimated by k1)
3. evaluate k3, the derivative at t¦¤t/2 (as estimated by k2)
4. evaluate k4, the derivative at t¦¤t (as estimated by k3)
All of these k values are samples of f(y,t) for likely-to-be-relevant conditions between t0 and ¦¤t. Accordingly, RK4 uses their (weighted) average to estimate an ‘average’ value of f(y,t) in that range.
RK4 works by estimating how the derivative f(y,t) changes on its way to the next moment, i.e. as yt goes to yt+¦¤t.
Let’s take a look at how it works. To do so, we walk through one iteration of RK4, where we are transitioning from t = 0, to t = ¦¤t. We will use a large ¦¤t=1, as picking a large ¦¤t exaggerates the visualization making things easier to see.
1. evaluate k1, the derivative at t0. 2. evaluate k2, the derivative at t¦¤t/2
(as estimated by k1)
3. evaluate k3, the derivative at t¦¤t/2 (as estimated by k2)
4. evaluate k4, the derivative at t¦¤t (as estimated by k3)
All of these k values are samples of f(y,t) for likely-to-be-relevant conditions between t0 and ¦¤t. Accordingly, RK4 uses their (weighted) average to estimate an ‘average’ value of f(y,t) in that range.
How well does it work?
Here we see RK4 more accurately integrating a function than Euler integration with a time-step 1/25th the size.
So even if RK4 took 25 times more calculations than Euler integration (it takes many fewer than that!), it still outperforms it.
Back to the Lotka- .
At the right we have a few different integrations of the lotka-volterra system.
Euler for different time step sizes and RK4.
Are the Euler trajectories “just a little bit off”?
Back to the Lotka- .
Errors are not just about being a little bit high or low in an estimate.
The systems we are interested in are non-linear, which means that a small difference can result in radically different dynamics.
An inaccurate numerical integrator can thus produce entirely different dynamics than those specified by the differential equation that it is trying to approximate.
Here is the same system simulated for a longer period of time.
The big-time step iteration diverges, with the amplitude of the oscillations growing and crashing more and more strongly, until we run into a numerical error / numerical artefact, i.e. a result of the numerics rather than the math they are trying to approximate.
Most obviously, we see negative population sizes for the prey around t = 17. These equations should never produce negative population sizes!
We can also see that the finer time step Euler integration amplitude is growing (more slowly), but growing nonetheless.0
This is indeed an ‘artefact’ of the numerics, i.e. not a real property of the differential equations, which specify a periodic oscillation (i.e. dynamics that precisely repeat every iteration).
CONCLUSION
Choosing your integrator and checking that your integrators are doing what they should be doing is very important.
Otherwise you might be studying a different system than you think you are studying!
RK4 is better than Euler integration, but it turns out it is still unsuitable for a set of dynamical systems termed stiff systems or stiff equations.
https://en.wikipedia.org/wiki/Stiff_equation
For these systems other methods such as ‘implicit’ integration methods must be used. These are more complicated and more computationally expensive.
Existing Libraries
Libraries exist in Python, matlab, C/C++, fortran and other languages for numerical integration.
library_solveivp.py
Existing Libraries
Libraries exist in Python, matlab, C/C++, fortran and other languages for numerical integration.
Why not just use them? Why did we have to learn all this crap?
Because if you don’t understand how numerical integrators work (at least to some basic level) it is hard to use them, or understand when they fail…
For example…
library_solveivp.py
For example if trying to simulate a bouncing ball, it might seem reasonable to write some INCORRECT code such as this for the derivative that describes the ball’s height…
But some integrators will hang/crash/fail when your differential equation(s) include discontinuities.
These change the step size during an integration to optimally compromise between the number of iterations and accuracy.
But they assume that the system is smooth. When they encounter a discontinuity, they try smaller and smaller step-sizes, in the assumption that eventually a ¦¤t that is small enough will make the system smooth.
Some integrators REQUIRE dynamics to be continuous, where as a bouncing ball’s dynamics involve discontinuities…
There exist ‘event’ mechanisms for dealing with discontinuities.
程序代写 CS代考 加微信: powcoder QQ: 1823890830 Email: powcoder@163.com