46

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


34

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 heuristic isn't ...


18

Pseudo time-stepping, probably better known as pseudo-transient continuation, is the technique of solving for the steady-state solution of time-evolving partial differential equations by setting an initial guess and using a time-stepper to evolve the solution forward. It tends to succeed where standard globalization strategies fail by taking advantage of ...


18

First thing, you could have mentioned, what RK method you have used. Here is a brief introduction to RK methods and Euler method, working, there merits and demerits. Euler method Euler's method is first order method. It is a straight-forward method that estimates the next point based on the rate of change at the current point and it is easy to code. It is ...


17

As Jed Brown mentioned, the connection between gradient descent in nonlinear optimization and time stepping of dynamical systems is rediscovered with some frequency (understandably, since it's a very satisfying connection to the mathematical mind since it links two seemingly different fields). However, it rarely turns out to be a useful connection, ...


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

While I haven't seen the exact formulation that you have written down here, I keep seeing talks in which people "rediscover" a connection to integrating some transient system, and proceed to write down an algorithm that is algebraically-equilavent to one form or another of an existing gradient descent or Newton-like method, and fail to cite anyone else. I ...


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

I do not think that there is a definite answer to this, because it might change from one topic to other (and also depends on the type of elements you are using). There are some recent papers talking about that, as well [2]. So, it is not a closed discussion. Furthermore, you can have different inertial components (at least in mechanics), when you have ...


12

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


11

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 is maintained quite accurately (I assume you mean to say that the orbits are elliptical, not circular). However, the method is still a numerical approximation ...


10

The reason is that with the exception of linear problems, if you do a Fourier (or other) decomposition in time, you end up with a significant number of problems that are coupled globally in time. In other words, you have to solve lots of problems on the entire time interval concurrently. That will typically bust your computational or memory budget. The ...


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

Already one good answer is available here, I just want to highlight some things, dual time stepping scheme uses pseudo time in addition to real time (so in your equations two time parameters will come one real and one pseudo). Real time act as a new dimension to the equations and generally discretised implicitly (refer paper). Real time step size should ...


9

A better way to look at it is that for a stiff problem, any stable explicit calculation leads to an error that is much smaller than the required error tolerance. There are many good methods for automatically detecting stiffness using explicit schemes, especially embedded Runge-Kutta pairs. See for example: Detecting stiffness with the Fehlberg (4, 5) ...


9

1. Can we numerically detect stiffness just by applying explicit methods? Suppose you have an initial value problem for some ODE on $[0,10]$. You take considerably large stepsize $\tau=1$ and an explicit Euler method, make your calculations with constant step size $\tau$ and get these points: You estimate the error and it appears to be big. Okay, then you ...


9

Following Kirill's suggestion, I ran the test with $N$ from a list of roughly geometrically increasing values, and for each $N$ computed the error as $$\epsilon(N) = \|\tilde z(2\pi) - \tilde z(0)\|_2 \quad\text{where}\quad \tilde z(t) = \left(\tilde{y}(t),\tilde y'(t)\right)$$ where $\tilde z$ represents the approximation obtained by numerical integration. ...


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


8

The PFASST (Parallel Full Approximation Scheme in Space and Time) and PEPC (Pretty Efficient Parallel Coulomb) algorithms have recently been used together to achieve parallelism in both space and time. PFASST does the time parallelism, PEPC does the space parallelism. The results of this were recently presented at the DD21 conference, and we have prepared ...


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

The Euler method does not take into account the curvature of the solution, so it tends to give different results depending on the step size. RK, depending on the order, takes into account the curvature. And this makes the estimated "next step" more accurate. Basically, if you are pretending a straight line is a good approximation of a curve (Euler) you will ...


7

I think it may be more useful to think of this is as numerical integration of a series of data points rather than as the solution of an ODE. Adams-Bashforth could work as suggested by @Omnomnomnom, but I think there are better methods. It seems that your acceleration values are "given" meaning that you have no control over the time step or the times at ...


7

It's a bit easier to see if you write your equation in the a semi-discretised system of the form $u^{\prime}(t) = F(u(t))$ and with the application of the $\theta$-method and approximating $u^{\prime}(t) \approx (w^{n+1} - w^{n})/\tau$ this gives, $$w^{n+1} - w^{n} - (1-\theta) \tau F(w^n) - \theta\tau F(w^{n+1}) = 0$$ with unknown vector $w^{n+1}$ and ...


7

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. If the integrator starts to use smaller step sizes then the process slows down. Looking at the particular form of your ODE, the $\sin(wt)$ factor has a ...


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


6

These methods can be roughly described in terms of two time-stepping methods, denoted here by $G$ and $F$. Both $G$ and $F$ propagate an initial value $U_n \approx u(t_n)$ by approximating the solution to $$ u(t) = u_0 + \int_0^t f(\tau,u(\tau)) \,d\tau $$ from $t_n$ to $t_{n+1}$ (that is, $\dot{u} = f(u,t)$). For the methods to be efficient, it must be ...


6

If the right hand side were independent of $u$ then one would generally use the averaged form $$ (1-\theta)s_(x,t^{n+1}) + \theta s(x,t^n). $$ In the nonlinear case you can't do that easily, as you note, but you can at least use some kind of extrapolation, for example approximate $$ (1-\theta)s(x,t^{n+1},u^{n+1}) + \theta s(x,t^n,u^n) \approx (1-\...


6

To answer your questions: As far as I know, in practice, if explicit methods require extraordinarily small time steps relative to your time scale of interest (see answers to this question on what it means for an ODE to be stiff) in order to yield accurate results, then for all intents and purposes, your problem is stiff. To determine requirements on step ...


6

In a certain sense, @Geoff Oxberry is correct in saying that stability and preservation of quadratic invariants are not directly related. For instance, there exist explicit methods that will preserve energy for your problem (and they are certainly not $A$-stable). However, in another sense there is a relation between the two. Well-posedness First, note ...


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