36

Yes, there aren't too many resources on this for some reason. For a very long time, the standard goto was "just use BDF methods". This mantra was set in stone for few historical reasons: for one Gear's code was the first widely available stiff solver, and for another the MATLAB suite didn't/doesn't include an implicit RK method. However, this ...


35

Since I just finished optimizing a lot of them in a software, DifferentialEquations.jl, I decided to just lay out a comparison of the main Order 4/5 methods. The Fehlberg method was left out because it's commonly known to be less efficient than the DP5 method. Backstories Dormand-Prince 4/5 The Dormand-Prince method was developed to be accurate as a 4/5 ...


21

van der Houwen's statement is correct, but it is not a statement about all fifth-order Runge-Kutta methods. The "Taylor polynomials" he is referring to are (as you seem to know) just the polynomials of degree $p$ that approximate $\exp(z)$ to order $p$: $$P_p(z) = \sum_{j=1}^p \frac{z^j}{j!}$$ For the fifth-order polynomial, it turns out that $|P_5(i\...


19

There are thousands of papers and hundreds of codes out there using Runge-Kutta methods of fifth order or higher. Note that the most commonly used explicit integrator in MATLAB is ODE45, which advances the solution using a 5th-order Runge-Kutta method. Examples of widely-used high-order Runge-Kutta methods The paper of Dormand & Prince giving a 5th-...


17

If you're doing celestial mechanics over long time scales, using a classical Runge-Kutta integrator will not preserve energy. In that case, using a symplectic integrator would probably be better. Boost.odeint also implements a 4th-order symplectic Runge-Kutta scheme that would work better for long time intervals. GSL does not implement any symplectic methods,...


17

There seems to be quite a bit of confusion about how to apply multi-step (e.g. Runge-Kutta) methods to 2nd or higher order ODEs or systems of ODEs. The process is very simple once you understand it, but perhaps not obvious without a good explanation. The following method is the one I find simplest. In your case, the differential equation you would like to ...


15

Bounds That is still true. In Butcher's book, page 196, it says the following: In a 1985 paper, Butcher showed that you need 11 stages to get order 8, and this is sharp. For order 10, Hairer derived a family of 17-stage methods, but it's not known if one can do better. The same information is given in Section II.5 of Hairer, Norsett, & Wanner vol. I. ...


13

TL;DR: It depends on what kind of accuracy you need. Energy conservation does not automatically equal accuracy. Suppose, you want to simulate the solar system, and you are using a solver that – to use an extreme example – just rotates the entire system by some angle every second. These solutions obviously conserve energy, but they are blatantly incorrect. ...


13

Let's take a stochastic differential equation: $$ X_t = f(t,X_t)dt + g(t,X_t)dW_t $$ Here's a few different arguments which lead to intuitive understandings of why the mathematics behind the higher order methods is necessary. I will be discussing in terms of strong order, which is the same as saying "for a given Brownian motion $W(t)$, how well does the ...


13

The problem you are encountering is likely not a consequence of your choice of algorithm, but in fact a consequence of the resulting dynamical system after applying time reversal. Per the definition of an attractor, all points in some neighborhood of the attractor will converge to the attractor under the flow of the dynamical system as $t\to\infty$. However, ...


12

As long as you're using standard double precision floating point arithmetic, very high order methods aren't needed to get a solution with high accuracy in a reasonable number of steps. In practice I find that the accuracy of the solution is normally limited to a relative error of 1.0e-16 by the double precision floating point representation rather than the ...


12

The Benchmark Setup In the Julia software DifferentialEquations.jl we implemented plenty of higher order methods, including the Feagin methods. You can see it in our list of methods, and then there are tons of others you can use as supplied tableaus. Because all of these methods are put together, you can easily benchmark between them. You can see the ...


10

Just to add to Brian Borcher's excellent answer, many real-life applications admit highly stiff ODEs or DAEs. Intuitively, these problems experience nonsmooth, abrupt changes over time, so are better modeled using low-order polynomials spread finely over short step-sizes, as opposed to high-order polynomials stretched over long step-sizes. Also, stability ...


10

The answer is quite simple. You are already comparing apples and oranges in the first equation. Garbage in, garbage out. The equation $y'=y$ if written properly is $$dy/dx=y.$$ Do you see it now? To correct it, simply write: $dy/dx=ay,$ where $a$ is a constant and in our example, $a=1$ in units of $1/x$.


10

This area has been fairly well researched, you may check e.g. Ketcheson's review of such methods: https://doi.org/10.1016/j.jcp.2009.11.006 which does contain some low-storage Runge-Kutta methods for fifth and sixth orders.


9

Both GNU Scientific Library (GSL) (C) and Boost Odeint (C++) feature 8th order Runge-Kutta methods. Both are opensource, and under linux and mac they should be directly available from the package manager. Under windows, it will probably be easier for you to use Boost rather than GSL. GSL is published under the GPL license, and Boost Odeint under the ...


9

256 equations is a relatively small number. All of the usual integrators, such as those included in Matlab, Maple or Mathematica should have no real problem with equations of this size and should be able to return answers in a fraction of the time it would take an algorithm you would implement yourself, because they use sophisticated explicit/implicit and ...


8

pressure is on the right track, but I will elaborate a bit. You are solving a system of coupled ODEs, however it appears you are viewing this as two systems of two equations. You need to combine them so that your system looks like: $$ \left(\begin{array}{c} \dot{x}_1 \\ \dot{x}_2 \\ \dot{v}_1 \\ \dot{v}_2 \end{array}\right) = \left(\begin{alignat}{1} ...


8

I am assuming that you are setting the error tolerance at 1e-12. You are correct that when an adaptive scheme accepts the current step size, it assumes the 5th order scheme was, for all intents and purposes, the "correct" answer. However this is only when it accepts the current step. If the difference between the 4th and 5th order steps are too large, it ...


8

You're asking how to produce dense output from your Runge-Kutta method. There are a number of ways to do this (see e.g. Hairer/Nørsett/Wanner). As noted in that reference, if you don't want to do more function evaluations aside from those already done by your fourth-order Runge-Kutta method, the best you can hope for is a third-order interpolant. This is ...


8

You can apply linear stability analysis. That is, for given $u=(x,v)$ compute the linearization $Df(u)$ of the right hand side if the equation is $u‘=f(u)$. The problem is stiff if those differ by orders of magnitude. At a glance, I would not expect this. You can determine a good step size by running the problem again with half the size. If the results are ...


8

These are all standard questions discussed in most books on ODE solvers. I would recommend Hairer & Wanner.


8

Notice that $\hat C_*=(1-r F(h))\hat C_n$, but the sign in front of $r$ is lost when you use $\hat C_*$ inside $\hat C_{n+1}$. I'm looking at p.149 of Numerical Solution of Time-Dependent Advection-Diffusion-Reaction equations by Hundsdorfer and Verwer on Google books. (I'm going to use their signs.) The stability region of Heun's method is $$|g(z)| = \...


7

The RK4 method implicitly constructs a degree 3 polynomial interpolant, using the data $f(x_i)$, $f(x_{i+1})$, $f'(x_i)$, and $f'(x_{i+1})$ in each interval. This interpolant can be constructed rather easily and efficiently using a linear combination of shifted Hermite basis functions in each interval.


7

@DavidKetcheson's answer hits the big points: you can always construct methods of high enough order using extrapolation, that's a very pessimistic bound and you can always do a whole lot better, all the good ones are derived by hand (with the help of some computer algebra tools), no lower bound is known, and the highest order methods are due to Feagin. Given ...


7

Another good starting point is the excellent paper of Kennedy, Carpenter, & Lewis (KCL2000): https://doi.org/10.1016/S0168-9274(99)00141-5 My own paper focuses more on the mechanics of low-storage implementation and what it implies for method construction, while KCL2000 is heavily focused on optimization and testing of methods. They give an extensive ...


7

The three best sources of Butcher tableaus are, according to me, the reference books Numerical Methods for Ordinary Differential Equations, 3rd ed., J.C. Butcher, Solving Ordinary Differential Equations I, 2nd ed., Hairer, Norsett and Wanner and Solving Ordinary Differential Equations II, 2nd ed., Hairer and Wanner. You will have to browse through the ...


7

The 14th order methods due to Feagin can be found in DifferentialEquations.jl. An example using them with 128-bit floating point arithmetic is as follows: using OrdinaryDiffEq, DoubleFloats function lorenz(du,u,p,t) du[1] = 10.0(u[2]-u[1]) du[2] = u[1]*(28.0-u[3]) - u[2] du[3] = u[1]*u[2] - (8/3)*u[3] end u0 = [1.0;0.0;0.0] tspan = (0.0,100.0) prob = ...


7

As @WolfgangBangerth already commented, this is usually referred to as a differential algebraic equation (DAE). These have their own challenges and there are special numerical methods for them. In general, for a system of DAEs you might have fewer differential equations than unknowns, so the algebraic constraint(s) serve to close the system and identify a ...


6

There are so many Runge Kutta methods, including Dormand-Prince 45 Cash-Karp 54 Fehlberge (sic) 78 Is there any comparison between them? Well, sure. Here are some traits to compare: Is the method implicit or explicit? (All of your examples are explicit RK methods.) What is the order of convergence? Are there any embedded error ...


Only top voted, non community-wiki answers of a minimum length are eligible