# Tag Info

63

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

40

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

19

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

17

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 . So, it is not a closed discussion. Furthermore, you can have different inertial components (at least in mechanics), when you have ...

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

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

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

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

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

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

Full space-time discretization of time-dependent partial differential equations is indeed a thing. If you use a structured mesh in time (in the sense that the time discretization does not depend on space) and appropriate choice of trial and test functions, you can fit several standard time-stepping methods (Crank-Nicolson, implicit Euler or some Runge-Kutta ...

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

6

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 property of the integration, it's not a stepwise property. Because of this, local error estimates and local changes of $\Delta t$ do not necessarily preserve the ...

5

Stability and preservation of invariants are unrelated. Hairer does a nice job of proving how various methods preserve different kinds of invariants in his book Geometric Integration. Chapter III is the most relevant to your question, where he demonstrates, for example, that RK methods and linear multistep methods preserve linear first invariants.

5

In the following, I am basically rephrasing p. 24 of "MCMC using Hamiltonian Dynamics" by Radford Neal. At least for leapfrog integration, one can analytically calculate the maximum time step for quadratic Hamiltonians. For such systems, a leapfrog step is a linear mapping $(q(t), p(t) \mapsto (q(t+\epsilon), p(t+\epsilon))$. Stability then depends on the ...

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

There are several equivalent ways of implementing Dirichlet boundary conditions with the finite element method. I'll give a brief overview, but you'll probably need more details, which can be found in the book Understanding and Implementing the Finite Element Method. It may be easier to analyze this for the Poisson problem $-\nabla^2u = f$ than for the heat ...

5

Forgive me for adding a me-too answer, but I just couldn't resist including this page from Press, et al: "Numerical Recipes in C":

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