Hot answers tagged

62

Let me start off with corrections. No, odeint doesn't have any symplectic integrators. No, symplectic integration doesn't mean conservation of energy. What does symplectic mean and when should you use it? First of all, what does symplectic mean? Symplectic means that the solution exists on a symplectic manifold. A symplectic manifold is a solution set which ...


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

So there is a ton to say about this, and we will actually be putting a paper out that tries to summarize it a bit, but let me narrow it down to something that can be put into a quick StackOverflow post. I will make one statement really early and keep repeating it: you cannot untangle the efficiency of a method from the efficiency of a software. The details ...


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

Although this post is now two years old, in case someone stumbles across it, let me give a brief update: Martin Gander recently wrote a nice review article, that gives a historical perspective on the field and discusses many different PINT methods: http://www.unige.ch/~gander/Preprints/50YearsTimeParallel.pdf There is now also a community website which lists ...


17

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. ...


16

Note that $\pi/2$ is represented in double precision format in a way that is not exactly equal to $\pi/2$. It's only accurate to about 15 digits. Thus you're starting every so slightly away from the equilibrium position. Since the equilibrium is unstable, it will eventually start moving.


15

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 ...


15

To complement Chris Rackauckas answer, to state some of the mathematical nonsense as well as some stuff you almost certainly know, a dynamical system is Hamiltonian if there is a description with coordinates $\mathbf{p}$ and $\mathbf{q}$ and a functional, $\mathcal{H(\mathbf{p},\mathbf{q})}$ such that $$\frac{d\mathbf{q}}{dt}=+\frac{\partial \mathcal{H}}{\...


15

I think the two main points have already been made by Brian and Ertxiem: your initial value is an unstable equilibrium and the fact that your numerical computations are never really exact provides the small perturbation that will make the instability kick in. To give a bit more detail how this plays out, consider your problem in the form of a general ...


14

DifferentialEquations.jl library is a library for a high level language (Julia) which has tools for automatically transforming the ODE system to an optimized version for parallel solution on GPUs. There are two forms of parallelism that can be employed: array-based parallelism for large ODE systems and parameter parallelism for parameter studies on ...


13

100 equations is not a particular large system. There are certainly many good integrators for this out there -- starting with Matlab's ode45 which should have no problems with a system of 100 equations. The challenge with ODEs is not typically the size, but the character. For example, is your system stiff? If so, you may want to look at CVODE. Do you need ...


13

This might seem extreme, but this can be analysed exactly. Take the system $$ \dot x_1 = x_2, \qquad \dot x_2=-x_1, \qquad x_1(0) = 1, \qquad x_2(0)=0. $$ Let $X=(x_1,x_2)$ be the state vector, $\delta t$ the time step, and $X^+$ the state vector for the next time step. Then the implicit Euler scheme is $$ X^+ = \delta t\left(\begin{array}{cc}0&1\\-1&...


13

The best approach is to use an ODE solver that is guaranteed to conserve the norm of the initial condition, i.e., for which $\|y_n\| = \|y_0\|$ for all $n\in\mathbb{N}$. Such solvers exist, and are called geometric integrators, since they preserve geometric properties of the exact solution (in this case, that energy is conserved, i.e., $\frac{d}{dt}\|y(t)\| =...


13

I'll write the equation short as $$\ddot x(t)+c\dot x(t)=a(t,x(t))$$ to separate the "easy" linear parts from the non-linear and forcing terms. On the first method The claimed order of the method is two, the implemented order is one. This is due to the implementation of the first derivative as a one-sided difference quotient. Changing that to the ...


12

At least one difference is that in a system of ODEs, all the equations are differential, e.g.: $$ \dot{x}=f(x,y)\\ \dot{y}=g(x,y) $$ whereas the definition of DAEs that I'm familiar with includes some non-differential (i.e. algebraic) equations in the set, e.g.: $$ \dot{x}=h(x,y)\\ y=l(x,y) $$ where $l$ is non-trival, and its solution can't be easily ...


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

odeint from the SciPy library defaults to the lsoda integrator described here. However, any simple description of asymptotic computation time is impossible. The reason is many fold. First, let me describe the algorithm. A common multistep algorithm for non-stiff equations are the Adams-Moulton methods. While these are implicit, the Adams-Bashforth methods ...


11

Is the shooting method the only general numerical method for solving BVP of nonlinear ODE(s)? No. Most other methods consist of three parts: Discretization. This may be done with finite differences, finite volumes, finite elements (Galerkin or collocation), spectral methods, and so forth. This reduces the problem from an infinite-dimensional one to a ...


11

And, it is my understanding that the 4 and the 5 are for the order of the global and local error, respectively. Your understanding is wrong. The local error of a Runge–Kutta method of order $n$ is proportional to $h^n$. What ode45 does is to estimate the solution (of one step) with two Runge–Kutta methods with local orders of 4 and 5, respectively (hence ...


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

There are the geometric integrations written by Ernst Hairer & co: E. Hairer, C. Lubich and G. Wanner (2002): Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential equations. Springer Series in Comput. Math., vol. 31 E. Hairer and M. Hairer (2002): GniCodes - Matlab programs for geometric numerical integration for ...


9

I would like to know if there are any practical methods for obtaining an approximation to $\mathbf{x}(t_f)$ (where $t_f \in \mathbb{R}$ is some given final time) which is provably correct to $N$ [sic] digits. That all depends on your opinion of the practicality of interval arithmetic. There are validated integrators available, such as the COSY code out of ...


9

What one usually wants in this situation is to preserve a discrete analog of time symmetry: namely, if the time discretization is applied to solve first forward and then backward in time, the initial condition is recovered. This is true if the method is invariant under the following substitutions: $$ \Delta t \to -\Delta t$$ $$ a^{n+j} \to a^{n-j}$$ (here ...


9

Any general-purpose ODE solver should be able to handle this linear coupled system of ODE very easily, for example: scipy.integrate.ode CVODE from the Sundials solver suite; it appears to have Python bindings here, and perhaps there are others. This kind of thing is typically discussed in any textbook on numerical methods. In general, computing the matrix ...


9

There are many kinds of RK methods which have extensions to exploit linearity. They all use some form of exponential or Lie Group idea (again exponential) to do so. Thus they generally do some form of integrating factor method and then apply the Runge-Kutta method to the IF transformed equation. If you want to see a listing of some of these, you can check ...


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

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} ...


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