Re: Octave-maintainers Digest, Vol 84, Issue 59

[c.]

This discussion is better suited for the help mailing list, so I am moving it 
there.

On 23 Mar 2013, at 15:14, address@hidden wrote:

> First, it is very likely that the equation you are solving is sensitive
> to floating point errors.  That is to say, a small numerical difference
> (maybe in the last bit of the mantissa) might later on be greatly
> exaggerated.  This is a well known phenomenon in solving ordinary
> differential equations (ODE), and it would be hard to say which, if any,
> is giving the correct answer - MATLAB or octave.  Another equation which
> illustrates this more profoundly is the famous Lorenz attractor
> http://en.wikipedia.org/wiki/Lorenz_system

> Secondly, your method of solving the (ODE) appears to be the so called
> Euler Method http://en.wikipedia.org/wiki/Euler_method.  It is well
> known to be rather inaccurate (hence your need for a very small time
> step).  Both octave and matlab have much more sophisticated ODE solvers,
> called lsode and ode45 respectively.  Or if you want to write your own
> solver, I would suggest the Runge-Kutta Method of order 4 with a fixed
> step size http://mathworld.wolfram.com/Runge-KuttaMethod.html.  It is
> actually fairly straightforward to program, and it will give far more
> accurate answers than your program.

This case is a long-time simulation of an Hamiltonian system. 
Rather than a standard high-order explicit solver like RK45 an explicit 
symplectic integrator is more appropriate for this problem.

I have an half-finished implementation of the velocity-verlet method which is a 
second
order symplectic integrator, here is an example of how to use to solve the OP's 
problem.
This integrator needs a bit more polishing and documentation, it will probably 
be included 
in the next release of odepkg.

c.

Orbit11.m
Description: Binary data

Orbit11b.m
Description: Binary data

velocity_verlet.m
Description: Binary data