Hot answers tagged

55

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

I'm not aware of any recent overview articles, but I am actively involved in the development of the PFASST algorithm so can share some thoughts. There are three broad classes of time-parallel techniques that I am aware of: across the method — independent stages of RK or extrapolation integrators can be evaluated in parallel; see also the RIDC (revisionist ...


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


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


19

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


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

As your matrix is independent of $u$ the result is a matrix exponential times the intial vector. The standard discussion of relevent method can be found from http://scholar.google.at by searching for ''Nineteen dubious ways''. For the scaling-and-squaring algorithm (the least dubious one), see also http://blogs.mathworks.com/cleve/2012/07/23/a-balancing-act-...


15

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


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


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

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


13

See David Stewart's new (2011) book on this topic, Dynamics with Inequalities: Impacts and Hard Constraints. Coulomb friction problems are mentioned several times in the analysis chapters. Chapter 8 is devoted to numerical methods for non-smooth ODEs and DAEs. It mostly advocates fully implicit Runge-Kutta methods with special treatment of nonsmoothness. ...


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

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

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


12

This is a very broad question and I am going to give you some things to think about (some are already included in your post, but they are repeated here for completeness). Scope of Problems You need to define the interface of how to specify problems. Are you going to allow parameters that can be fixed or can vary for solutions? Are you going to allow ...


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

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


11

To get something that looks realistic for planetary orbits, you shouldn't use the forward or backward Euler methods. These will cause your planets to spiral outward or inward. You should use a symplectic method. You may also need to adjust the timestep to be smaller when two bodies are very close to each other. Read Chambers (1999) A hybrid symplectic ...


11

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)\| =...


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


11

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


10

The most significant reference I know of is David Stewart's thesis, which is more than 20 years old: High Accuracy Numerical Methods for Ordinary Differential Equations with Discontinuous Right-hand Side The abstract references several significant earlier works. A keyword here is differential inclusion.


10

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


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

You won't find one reference systematically covering the analysis of all the important methods for PDE. The field of discretization techniques for PDE is at least an order of magnitude larger than either topic you mentioned above. For any methods involving implicit solves, studying discretizations without also considering solution methods (e.g., associated ...


9

You seem to think that the matlab integrator is the problem and by this hypothesis you should try a different integrator. But before you do this I think it would be useful to check whether that's actually the case. For example, it may be that your user-provided function hangs. Or that it is just very slow. Or that it produces wrong results. As programmers we ...


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