83 votes
Accepted

What does "symplectic" mean in reference to numerical integrators, and does SciPy's odeint use them?

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 ...
Chris Rackauckas's user avatar
43 votes
Accepted

BDF vs implicit Runge Kutta time stepping

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 ...
Chris Rackauckas's user avatar
17 votes

What does "symplectic" mean in reference to numerical integrators, and does SciPy's odeint use them?

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 ...
origimbo's user avatar
  • 2,239
15 votes
Accepted

Why is leapfrog integration symplectic and RK4 not, if the latter is more accurate?

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 ...
Wrzlprmft's user avatar
  • 2,032
13 votes
Accepted

Why do I get an oscillatory solution when applying the implicit trapezoidal method to the linear diffusion equation?

I think you might have mixed up some terminology. An exponential integrator would use some type of eigensolver or related approach to exactly calculate $$f^{n + 1} = e^{\delta t\cdot \mathcal A}f^n$$ ...
Daniel Shapero's user avatar
11 votes

Why do planets move at the wrong speed in my solar system model?

This is expected behavior with the Verlet algorithm. It is a symplectic integrator, which means that it will preserve quadratic invariants to within roundoff error -- thus the form planetary orbits ...
David Ketcheson's user avatar
10 votes

Do there exist low-storage Runge–Kutta methods with an order larger than four?

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 ...
Reid.Atcheson's user avatar
10 votes
Accepted

Which numerical methods preserve time reversal symmetry?

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 ...
David Ketcheson's user avatar
10 votes
Accepted

Why is velocity Verlet better than Verlet for gravity if it has a worse order of magnitude for the error term

That is a common misconception. Verlet in whatever form is a second-order method. The misconception results from the fact that the truncation error of the Verlet formula is of 4th order. The naive ...
Lutz Lehmann's user avatar
  • 5,984
8 votes
Accepted

Computational time not proportional to integration interval in ODE-solver?

Substantially edited, since the original poster changed his equation... In general, the MATLAB (and Octave) ODE solvers dynamically adjust the step size as needed to maintain an accurate solution. ...
Brian Borchers's user avatar
7 votes

Do there exist low-storage Runge–Kutta methods with an order larger than four?

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-...
David Ketcheson's user avatar
7 votes
Accepted

The velocity Verlet method and variable time steps

The second formula is just the velocity verlet, and it's correct but if you adapt time steps then it's not symplectic. In a separate answer I describe in quite detail that symplecticness is a global ...
Chris Rackauckas's user avatar
7 votes
Accepted

Solving the time dependent Schrödinger equation with leapfrog integration in 1D

I read through the paper you linked and they give the stability condition for this method to be (eq. A6) $$ \frac{-2}{\Delta t} \le V \le \frac{2}{\Delta t} - \frac{2}{m \Delta r^2} $$ This has to be ...
helloworld922's user avatar
7 votes
Accepted

Does this second-order implicit Runge-Kutta method have a name?

As Wolfgang stated in the comments, this is not a traditional RK due to the inconsistent time evaluations within a stage. At first it would seem it can't even be cast as an additive RK since terms ...
Steven Roberts's user avatar
6 votes

Comparison of velocity Verlet and leapfrog algorithms

The shortest answer is that the Verlet / Leapfrog methods are symplectic and time reversible, and these are desirable numerical properties that reflect the physical reality of certain simulation ...
Richard Zhang's user avatar
6 votes
Accepted

Solve an ODE with positivity-preserving property unconditionally

The two properties are usually called positivity-preserving and monotonicity-preserving (makes it easier to find this question). Looking at http://www.ams.org/journals/mcom/2006-75-254/S0025-5718-05-...
Kirill's user avatar
  • 11.4k
6 votes
Accepted

How does non-dimensionalization improve the behavior of ODE solvers?

Background Let’s have a brief look at what parameters an integrator typically has and what they are used for. All of them govern step-size adaption: Relative tolerance (...
Wrzlprmft's user avatar
  • 2,032
5 votes

Test of 3rd-order vs 4th-order symplectic integrator with strange result

Plotting the error of $\ddot q=-q$, $q(0)=1,\dot q(0)=0$ over the full interval, scaled by the power of the step size given by the expected order, gives the plots As expected, the graphs for ...
Lutz Lehmann's user avatar
  • 5,984
5 votes
Accepted

Fourth order IMEX Runge-Kutta method

I think the work of Kennedy and Carpenter (mentioned already by @GoHokies) is still the definitive study on this topic. The journal paper can be found here; for some reason Google Scholar only ...
David Ketcheson's user avatar
5 votes
Accepted

Developping PDE with Python symbolically and numericaly

But how do you numerically (Runge-Kutta or whatever) solve the objects coming out of SymPy? That depends a lot on your specific problem and what level of optimization you need. Most ODE solvers (...
Wrzlprmft's user avatar
  • 2,032
5 votes

Time integration of wave equation

It's not really meaningful to talk about integrating the equation in form A or B, since one way to integrate A is to first transform to B and then discretize. You can only really compare the actual ...
David Ketcheson's user avatar
5 votes

Numerical integrator for $a'(t)=e^{-a(t)}f(t)$

You want a numerical solution, but this might help you check your computed results. If $a$ satisfies the ODE, you know $e^{a(t)}a'(t) = f(t)$. Integrating you get \begin{align} \int_0^t\, f(\tau)\, d\...
A rural reader's user avatar
4 votes

How do you apply boundary conditions in a time-stepping problem?

As HBR mentioned, the boundary conditions can often be immediately incorporated into $A$ and $b$. For example, suppose we wish to solve the 1D heat equation with Dirchilet boundary conditions $$ u_t =...
eepperly16's user avatar
4 votes
Accepted

Symplectic Algorithms for Hamilton’s Equations as opposed to just Volume-Preserving

Backward error analysis comes to mind. For example, for the Verlet scheme, you can say that the numerical solution turns out to be the exact solution for a perturbed Hamiltonian system (also known as ...
Juan M. Bello-Rivas's user avatar
4 votes
Accepted

How to select initial time step in adaptive time step ODE solver (TR-BDF2)

Many numerical tips and theoretical explanations can be found in this book from Hairer and Wanner: https://www.springer.com/gp/book/9783540566700 In this book, a strategy is described, which uses a ...
Laurent90's user avatar
  • 1,808
4 votes
Accepted

Write incompressible Navier Stokes as ODE in $(\mathbf{u},p)$

One approach to convert this into an ODE is with index reduction methods. These allow you to convert high-index DAEs into low-index DAEs or ODEs. See section VII.2 of "Solving Ordinary ...
Steven Roberts's user avatar
4 votes
Accepted

How to deal with solving coupled ODE systems where variables are updated multiple times within each timestep?

It would have been helpful to see the ODE you're actually trying to solve, but the way I interpret what you write is that the right hand side of your ODE system consists of a number of terms ...
Wolfgang Bangerth's user avatar
4 votes

Using velocity verlet algorithm for nbody simulation results in planet leaving orbit

There are two sources for such an error, if the exact solution is known to stay bounded. Verlet becomes catastrophically incorrect in singular situations, that is, if two objects become close in a N-...
Lutz Lehmann's user avatar
  • 5,984
3 votes

Developping PDE with Python symbolically and numericaly

I want to point out that PyODESys is a Python library that does exactly what you're looking for: send SymPy expressions to ODE solvers. It's nice because it links to more solvers, but @Wrzlprmft's ...
Chris Rackauckas's user avatar
3 votes
Accepted

Open source solver for continuous-time stochastic non-linear DAEs (SDAEs)

DifferentialEquations.jl in Julia can do it if you can write it in mass-matrix form. You won't find it mentioned in the tutorial, but you can provide a mass matrix as part of the ...
Chris Rackauckas's user avatar

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