Re: [Phys-L] accurate numerical solution of equations of motion

[John Denker <jsd@av8n.com>]

On 04/03/2014 09:38 PM, Bernard Cleyet wrote:
> I’ve been using the leapfrog method described in Eisberg and Lerner,
> I thought successfully.  Until I tried testing the accuracy of a very
> simple harmonic oscillator by fitting and then having the fit find
> the residuals.  The residuals increase with time!   A quick read of
> Cromer explains and suggests I switch to the last point method.
> However, my understanding of this conversation is I will have error
> when I introduce drive and dissipation.   I’ve already run these, but
> only for a few seconds of a one second pendulum.  Please assist.
>
> Note: the error(s) I found after 80 periods of my one second harmonic
> oscillator with an amplitude of 0.1 radian are ~ + & - 2e-6 radian
> [time step = 2E-3s;  difference between a sinusoidal fit and the
> simulation]

Leapfrog is a second-order method.  It is a /symplectic/ 
method, which means its long-term behavior is extraordinarily
good for a /simple/ harmonic oscillator or other Hamiltonian
system.  The errors should be dominated by machine precision
(one part in 10^15 per step) rather than by series truncation
error in the equations of motion (stepsize squared).

Symplectic structure does not however guarantee that 
the  frequency is correct.  Instead it guarantees that
/at every zero crossing/ of the numerical model the 
energy is correct, to very high precision.  If the
numerical model is out of phase with the exact solution,
you need to change the numerical method somehow.

There exist high-order symplectic methods.  Based on the
limited available information, my guess is that fourth order 
should be plenty for this application.

Canned routines for doing this exist.  They are not available
in True Basic, but they are available.
   http://www.av8n.com/computer/geometric-numerical-integration/
   http://www.av8n.com/computer/geometric-numerical-integration/README

I recommend sticking with symplectic methods.  A fourth-order
symplectic method won't be any worse than a non-symplectic 
method of the same order.  For a Hamiltonian system it will be
much better.  For a /slightly/ dissipative system it should
retain some of its advantages.

The conceptual basis for symplectic integrators is discussed 
here:
  http://www.av8n.com/physics/symplectic-integrator.htm

It is profoundly related to the physics of Liouville's theorem
  http://www.av8n.com/physics/liouville-intro.htm
which is in turn related to the second law of thermodynamics,
the Heisenberg uncertainty principle, the brightness theorem 
(in optics), the fluctuation/dissipation theorem, et cetera.

===============================

I carry a pocket tool that includes a tiny pair of pliers.
I also have a lug wrench in the trunk of my car.  I use the
tiny pliers 100 times more often than I use the lug wrench.

However, when it comes time to tighten lug nuts, I use the
lug wrench.  If you try to do this using pliers, you'll 
just hurt yourself and strip the nuts long before you get 
them tight.

I mention this because integrating the equations of motion
is a lug-nut type problem, and True Basic is a pocket pliers
sort of tool.  This is the smallest and easiest lug-nut you
are ever going to see, but it is still outside the reach of
the tiny pliers.

Learn how to use the C++ compiler.  Trying to do without
is penny wise and pound foolish.  If you're stuck with a 
Windows box, install cygwin.
  http://www.cygwin.com/

C++ is designed to be extensible.  Basic comes with sine
and cosine built in, but if you want to use some support
routines that aren't built in, things get real complicated
real fast.  C++ has more of a barrier to entry in the short
run, but it is easier and better in the long run and even
in the medium run.  There are huuuge numbers of support
routines that you can use, without having to type them in
yourself.
  Hint:  boost
  Hint:  gsl
  Plus lots of other stuff.