Often different numerical algorithms are judged based on how easy they are to implement, order of error and number of functional evaluations required. In "most cases" the computational speed is determined by the step size, and number of functional evaluations required. And, for most "regular" functions RK-4th order has been well acknowledged for its efficiency and suitability. May be that is why in many software packages RK is used as the default ode solver.
Can anyone suggest a numerical algorithm which can give a more accurate result than RK-4th order? I am ready to compromise a little (say 10-20% more), on the computational speed (as it can be circumvented using parallel programming, multi-threading or GPU stuff).
Another query:
Also, numerical integration algorithms which are better than the typical Simpson one-third or Gaussian 4-point quadrature rules?
There are several methods rather than Runge-Kutta for example "Adams-Bashforth-Moulton", see, https://en.wikipedia.org/wiki/Linear_multistep_method
"Euler method", see, https://en.wikipedia.org/wiki/Euler_method
and others, each one can solve the ODEs according to the stiffness of the ODE, the Runge-Kutta, for example, can solve non stiff ODEs while the Trapezoidal rule can solve stiff ODEs.
If you are working with Matlab then there is a lot of solves like, ode23t, ode15, ... you can solve your problem by changing the solver until you reach the rigth solver (i.e. the behavior of the results, and the time).
Best regards
Muhammad
@abdulrasool
No idea why you suggest Euler and trapezoidal method. Are you aware that order of error in both cases is huge enough, and so the solution obtained is often considered unacceptable?
There's a huge variety of methods out there and the type one selects depends on variety of factors such as whether the system is stiff or not. I think it is true that RK4 is the most accurate RK scheme which is of n^th order with n intermediate steps. A fifth order RK method needs six steps, for example. With the computations I do for ODEs, I rarely need to go beyond a 4th order RK scheme. It might be worth doing some numerical experiments to compare the work done by a selection of methods in order to attain the same absolute accuracy. So a 4th order method requires fewer timesteps than a 2nd order method to attain a predetermined accuracy, but each timestep uses more algebra. I don't know the answer, but I suspect that it is easy to do.
With regard to numerical integration, why not turn it into an ODE? If you need the integral of f(t) between x=a and x=b, then solve dy/dt=f(t) subject to y(a)=0, and then your integral is y(b).
But if you want a nice explicit formula, then Boole's rule has 6th order accuracy. It is possible to derive it by applying Richardson's extrapolation to Simpson's rule (and also derive Simpson's rule from the Trapezium rule). I set this as an exercise for my year 1 undergraduates this year for the first time!
Vikash,
Back in the early 80's when I was a graduate student, I took a course in the numerical solution of ODE's. One of the assignments was to compare a variety of methods for a suite of (three) benchmark ODE's. There were six algorithms we compared, two of which were fourth order RK methods (one that existed as a canned routine and another that we programmed ourselves). The IVP's were of the form y' = f(t,y) for t in [0,b]. The three example ODE's were f(t,y) = 2ty, Y0=1; f(t,y) = -2ty, Y0=1; and f(t,y) = t^3.5, Y0 = 0. The values of b were 1, 4, 7, and 10 and step sizes were 1/4, 1/16, 1/64, and 1/256.
One of the better methods (for accuracy), if not the best, was a predictor/corrector method where the predictor is a fourth order Adams-Bashforth algorithm and the corrector is a fifth order Adams-Moulton algorithm. The corrector was applied once for one of the algorithms and twice for another one. Unfortunately, I can't provide you with timing results, but I'm guessing that timing results from an 1980's machine may not be very relevant today anyway.
This is of course just one data point and the more general remarks that Prof Ree's has already made I'm in complete agreement with,
@ Wilson
I will try programming the method you have suggested. Good info though!
Considering that, you need computational speed, so I think that, an adaptive method can be applied. For instance, MATLAB function ode45 can be used. This function is based on an explicit Runge-Kutta (4, 5) formula with a variable time step for efficient numerical computation of differential equations.
The predictor-corrector methods are, generally, faster than corresponding-order Runge-Kutta methods, e.g. 4th order predictor-corrector methods are nearly twice as fast as 4th Runge-Kutta methods. Both methods are comparable in accuracy.
Dear Fabrizio, which P-C method do you have in mind? A standard fourth order Runge-Kutta scheme uses four function evaluations per timestep, whereas, for example, the fourth order Adams-Moulton scheme (as given by equation (3.34) of http://www.cems.uvm.edu/~tlakoba/math337/notes_3.pdf ) has five function evaluations per timestep. Superficially this looks as though the P-C method is slower than the RK method if the timestep is the same. But the P-C method is implicit and may be used for timesteps for which the RK method is unstable, so it could be faster. So my view is that it is difficult to compare these different methods in general, and much depends on the equation being solved and how stiff it is
Andrew
For example, assume a fixed step size, then the 4th order Adams-Moulton predictor-corrector (without repeated corrector iterations) requires half THE derivative evaluations as THE 4th order Runge-Kutta methods. In this case the predictor-corrector methods are "presumed" to be 2 times faster than RK.
Moreover, PC methods with variable steps are not unconditionally faster than RK methods.
Fabrizio, I now understand. If we are solving y'=f(y,t), and if one is at timestep t_i and wish to calculate y_{i+1} , then the Adams-Moulton P-C method has already computed values of f_{i-1}, f_{i_2} and f_{i-3} in previous timesteps. So in order to find y_{i+1} one will have to evaluate f_i and the predicted value of f_{i+1}. So it is two new RHS evaluations, as opposed to four for RK4. Nice. Perhaps I need to start using the Adams-Moulton method from now on! A happy day! Thanks for making me think about this.
Thank you prof Andrew S.Rees for the informations , he can use also the RK of the ordre of 5
I had a quick look at the predictor-corrector (PC) methods. I am wondering about your claims that it could be faster than RK4. Though, it is true that the order of error in RK4 is the same as that of PC.
For PC to work, there are many more functional evaluation and algebraic calculation embedded in predictor part as well as in the corrector part. When combined together, it has to go many more steps than RK4. Then how can it be faster?
I will code it and check.
Re Vikash's post on July 9th regarding the relative number of function evaluations for PC methods as opposed to RK methods. At first glance they look as though RK4 has four per timestep according to the standard formula, while the usual 4th order PC method has five. However, three of those five have already been evaluated during earlier timesteps, and therefore only two new ones need to be found. On the other hand, one must evaluate all four for the classical RK4 method. This is why I said, "Happy days", in my previous post(!), and it is also why I need to re-evaluate the near-blanket use of RK4 as my computational workhorse for ODEs over the last few decades.
Andrew
Gear's method, implemented in Matlab as ode15s and in SciPy as method='bdf' , is better (more stable) on stiff systems and faster on lower order systems than Runge Kutta 4-5. On higher order non-stiff systems, it is about the same as RK45, implemented in matlab as ode45 and in NumPy as ode45.
I tried Google to find the algorithm for Gear's method for different order (if there is?), but did not found any satisfactory result. Can anyone send or provide a link for the same?
Past couple of months I have been doing exhaustive testing of numerical algorithms on different sets of differential equations and this is what I have found. I have tested algorithms from implicit families: Euler to Adams-Moulton, and also explicit families: Euler to RK8 Order. Please feel free to comment on them:
1) In general, higher order explicit RK family algorithms give good solutions for most problems. In some sense, they can be considered reliable too. But NOT completely, at least for me.
2) The higher the RK order, the lower is the error, is NOT always true. Example, RK6 order will not always give better solution than RK4 or even RK3. That may differ from case to case. Though in general, higher order RK gives better solution than the lower order ones.
3) But there are so called stiff systems, for which RK families could not give stable solution unless the grid-spacing is very very fine, i.e., extreme computational load becomes necessary. I have encountered problems in which I had to fine-grid x from 0 to 2 in 2500 points for RK4 to give decent solution.
4) In contrast, for stiff problems, implicit families of Euler, Adams-Bashforth and Adams-Moulton often yield stable and accurate solutions with much coarser grid. That is, I achieved RK4 type results with just 1000 points instead of 2500 using Implicit Euler.
5) I have not tested the implicit algorithms for those problems on which explicit algorithms give good solutions at ease. So, I have to find drawbacks of implicit algorithms, if any.
6) It is often difficult to guess by the appearance of an equation if they are stiff or not. I had an initial idea that presence of exponential or trigonometric terms lead to stiffness but this is not often true.
7) As I said above for stiff problems, explicit families give stable solution with very fine grid. But I have also noticed that with more refinement of grids, the "nature" of the solution may not change, but there is always some minor shift in the solutions observed.
8) So, in conclusion, in order to make sure that the solution obtained is good enough, the problem should be solved by at least two algorithms from implicit and explicit family each and check if they yield similar nature of solutions, if not exactly identical. This is crucial specially when the true solution is not available for comparison, and yet the results are going to be a key part in a prospective publication.
There is a new method called piecewise analytic method (PAM).
Preprint Piecewise Analytic Method VS Runge-Kutta Method (Comparative Study)
Piecewise analytic method (PAM) is a new technique used for solving nonlinear differential equation. The advantages of PAM are( these points are proofed):
1. It gives a general analytic form that can be used in differentiation and integration.
2. It can solve highly non-linear differential equation.
3. The accuracy and error can be controlled according to our needs very easily.
4. It can solve problems which other famous techniques can’t solve.
5. In some cases, it gives the exact solution.
Preprint Solving Nonlinear 2 nd Order Differential Equations Using Pi...
Preprint Piecewise Analytic Method (PAM) is a New Step in the Evoluti...
Dear Vikash and All,
If you are interested in fast ODE solver on GPUs, plesase visit my website: www.gpuode.com.
I have a Research Gate project as well:
https://www.researchgate.net/project/Integration-of-independent-ODE-system-on-GPUs
All the best,
Ferenc
Dear All,
Is there a method can compete Runge Kutta method?
You can see the answers :
https://www.researchgate.net/post/Is_there_a_method_can_compete_Runge_Kutta_method
Dear All, I have already solved a system of strongly coupled nonlinear ODEs (Keller-Miksis equations) using RK algorithm in Matlab. Here it is
https://journals.aps.org/pre/abstract/10.1103/PhysRevE.99.042209
But I have also felt that for more complex multi-body problems in larger computation domain, we better exploit particularization of CPUs (most workstations comes with 8 cores) or use an even faster GPUs.
Also, RK method is really well established and reliable for initial value problems that are encountered in cavitation bubble dynamics. Adams-Bashforth and Adams-Moulton are much better algorithms for stiff problems though.
Another thing that I noticed is that Matlab ode solver uses RK algorithm plus something "extra" in its solver algorithms - some sort of conditioning I guess. For some sensitive problems, RK when hard coded in Fortran or Python could not give any result. I guess this was because of emergence of some imaginary numbers in its intermediate calculations.
A good approach that I believe is that solve the same problem using explicit algorithms such as RK or its adaptive variations, and implicit algorithms such as Adms-Moulton like algorithms just to make sure the"nature" of obtained solution is at least consistent.
Dear All,
I agree with Vikash, and thanks for sharing your paper. It is really interesting for me.
I only have some minor comments. The main strength of my GPU implementation is actually the framework with which one can exploit the high processing power of GPUs. We already working to extend it with other numerical schemes and with efficient coupling of identical systems.
Another point is, if one starts to use massively parallel hardware architecture (e.g. GPUs), the choice of best algorithm for your problem becomes highly non-trivial. Even if you have a stiff problem, non-stiff solvers can outperform the stiff ones. The reason is a bit complicated and it can depend on problem to problem, it is the consequence of how the hardware handle thousands or millions of threads. Have a look on the documentation (Performance consideration) of my code in its github repository:
https://github.com/FerencHegedus/Massively-Parallel-GPU-ODE-Solver/tree/master/PerThread
My articles and ideas:
Article A New Reliable Numerical Method for Computing Chaotic Soluti...
Article A New Accurate Numerical Method of Approximation of Chaotic ...
Preprint An Accurate Numerical Method and Algorithm for Constructing ...
Method Solution of the Duffing Equation by the Power Series Method ...
By the quality of approximation, this modification of the power series method is very good.
i'v Java Application for solving differential equations using Runge-Kutta fourth order , for more : https://www.researchgate.net/publication/340871252
Almost Runge-Kutta methods (ABDULRAHMAN NDANUSA & ZAINAB ISHAQ ABDULAZEEZ at department of mathematics, federal university of technology minna, niger State, Nigeria) email: [email protected]
1. RK4 is the simple fixed-step 4th order method developed by W. Kutta in 1901 following the pattern of K. Heun's methods from the previous year [1]. Even in the original paper it was observed that, while it leads in simplicity, the error of the 3/8 method is smaller by about a half. [2]
- [1] https://math.stackexchange.com/questions/2527302/whats-the-motivation-for-runge-kutta-methods/2527316#2527316
- [2] https://math.stackexchange.com/questions/3167436/why-use-classic-fourth-order-runge-kutta-over-the-3-8-rule
2. The algorithm that you refer to as default solver implemented in most solver packages (e.g., as ode45, dopri5 or RK45) is the Dormand-Prince embedded (4)5 (FSAL - first same as last) method (1980) [3]. It is a 5th order method with automatically locally adapted step size. It is still simple to implement [4] if you do not want to use the library version, one only has to be careful to correctly translate the coefficients of the Butcher tableau. For an overview of frequently used (4)5 methods see [5]
- [3]Article A Family of Embedded Runge-Kutta Formulae
- [4] https://stackoverflow.com/questions/54494770/how-to-set-fixed-step-size-with-scipy-integrate/54502790#54502790
- [5] https://scicomp.stackexchange.com/questions/24998/which-runge-kutta-method-is-more-accurate-dormand-prince-or-cash-karp/25031#25031
3. What is the most efficient method depends also on the problem to be solved, but as far as I understand the experts [6], the optimum is somewhere between 7th and 10th order (embedded) methods, in some benchmarks even 12th order methods are competitive.
- [6] https://scicomp.stackexchange.com/questions/25581/why-are-higher-order-runge-kutta-methods-not-used-more-often/25619#25619
- Badges
- Science topic