The Stack Overflow podcast is back! Listen to an interview with our new CEO.

Hot answers tagged

45

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


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


32

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


27

There is a wide variety of algorithms; Barnes Hut is a popular $\mathcal{O}(N \log N)$ method, and the Fast Multipole Method is a much more sophisticated $\mathcal{O}(N)$ alternative. Both methods make use of a tree data structure where nodes essentially only interact with their nearest neighbors at each level of the tree; you can think of splitting the ...


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


17

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


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

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

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

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


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

Look at the fast multipole method. It is highly scalable and $O(n)$. It allows trading off between precision and cost. Here's an example where it is run at 42 Tflops on a GPU cluster.


13

It looks like the equations you're dealing with are all polynomial after clearing denominators. That's a good thing (transcendental functions are often a little harder to deal with algebraically). However, it's not a guarantee that your equations have a closed-form solution. This is an essential point that many people don't really "get", even if they know it ...


13

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


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

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


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


10

Stiffness involves some separation of scales. In general, if you are interested in phase of the fastest mode in the system, then you have to resolve it and the system is not stiff. But frequently, you are interested in the long-term dynamics of a "slow manifold" rather than the precise rate at which a solution off the slow manifold approaches it. Chemical ...


10

Part 1 Small eigenvalues are not included in the definition of stiffness for ODE (initial value problem) systems. There is no satisfying definition of stiffness that I know of, but the best definitions I've come across are: If a numerical method with a finite region of absolute stability, applied to a system with any initial conditions, is forced to use ...


10

You might consider it "overkill", but PETSc's time integration package can be used with C99 complex (configure --with-scalar-type=complex). Supported methods include explicit Runge-Kutta low-memory strong stability-preserving Runge-Kutta Rosenbrock-W additive Runge-Kutta IMEX These implementations are most appropriate for high-dimensional problems such as ...


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

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

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


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