# Tag Info

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

37

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

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

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

13

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

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

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

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

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

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

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

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

6

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-01794-1/S0025-5718-05-01794-1.pdf and http://homepages.cwi.nl/~willem/DOCART/SIAM_HRS.pdf it seems that implicit Euler is the exception among implicit linear ...

5

I believe you're switching around the definitions of "inner" and "outer" stepping (this nomenclature gets worse in segregated schemes where you have at least 3 separate iterative schemes going on.) The core of what you're missing is that "pseudo-time" isn't a time at all; it's just a convenient way to achieve an iterative solution. Take a 1st order backward ...

5

The following snippet of code is an implementation of RATTLE on a system with the constraint $g(x, y) = K x^2 + y^2 - 1 = 0$. // Constraint. double g(const double2& r) { return K * r.x * r.x + r.y * r.y - 1.0; } // Gradient of constraint. double2 G(const double2& r) { return double2(2.0 * K * r.x, 2.0 * r.y); } void rattle(double2& q, ...

5

For something with a spectral flavor in time, look at deferred correction methods, starting with this paper. I would argue that they're not spectral in the usual sense of the word, but they give you a family of arbitrary-order Runge-Kutta methods, so if you think of "refining" by increasing the order (by adding more nodes), then the convergence can be ...

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