Re: [Phys-L] pendulum +- numerical methods

[John Denker <jsd@av8n.com>]

Keep in mind that if all you want to know is how the
nonlinearity affects the pendulum, you can express that in
terms of elliptic integrals.  For not-too-large angles the
function is easy to evaluate.  Any scientific math library 
worthy of the name will have a canned function to do it.
See also figure 3 and the formulas at
  
https://en.wikipedia.org/wiki/Pendulum_%28mathematics%29#Arbitrary-amplitude_period
or just
  https://www.google.com/search?q=%22elliptic+integral%22

On 09/23/2014 03:02 PM, Bill Nettles wrote:
> I've taken a look at the code, and I don't see how you are
> calculating what you call "residuals."  Is that a root mean square
> error of a fitting routine, the root mean square of the
> point-by-point comparison to a cosine, or the max difference of any
> one calculated point compared to a cosine, or something else?
> 
> If I take sqrt[ 1/N *sum( (calc-cosine)^2) ] I get numbers in the
> range of 2.7e-7 for 0.1 ms time step over 50 periods.  I get 1e-5 for
> 1.0 ms time steps over 50 periods.

We agree that's a weird thing to calculate in this situation.
The mean-square residual over 50 periods is heavily dominated 
by the phase error in the later cycles ... and in that case 
you'd be much better off measuring the phase error directly.

Actually it scales like the square of the phase error.

Using the second-order symplectic approach i.e. Velocity Verlet,
  https://www.av8n.com/physics/symplectic-integrator.htm#sec-verlet
I get the following values at the end of 50 cycles:

  t/period          x              v          Delta E      rms resid       ms 
resid      steps/cycle
 50.000000   9.986638e-01  -5.165156e-02  -2.635697e-06   2.110231e-02   
4.453074e-04            100
 50.000000   9.999999e-01  -5.167710e-04  -2.620792e-12   2.109700e-04   
4.450836e-08           1000
 50.000000   1.000000e+00  -5.167713e-06   3.086420e-14   2.109840e-06   
4.451427e-12          10000
 50.000000   1.000000e+00  -5.167711e-08  -1.345590e-13   1.057928e-08   
1.119212e-16         100000
 50.000000   1.000000e+00  -5.168329e-10   7.371881e-14   1.023437e-07   
1.047424e-14        1000000

Note that the "v" value is also the phase error (in radians).

The second-order method is orders of magnitude more accurate
than any first-order method ... yet is super-simple to implement.

Note that in the last two rows, the "residuals" are so tiny 
that they are getting clobbered by machine roundoff errors,
but the "v" values are just fine, and continue to scale
like the /square/ of the step-size, as we would expect for
a second-order algorithm.

In all rows except the first two, the energy error is so small 
that it is getting clobbered by machine roundoff.  It appears
to scale like the sixth power of the step size.  It's not
obvious to me why that is;  I would have guessed fourth power.

Also note that you don't have to grind it out to 50 cycles.
The phase error is pretty small after the first 1 or 2 cycles:
  t/period          x              v          Delta E      rms resid       ms 
resid      steps/cycle
  1.000000   1.000000e+00  -1.030812e-11  -1.283418e-13   1.355260e-11   
1.836731e-22        1000000
  2.000000   1.000000e+00  -2.067598e-11  -3.863576e-14   2.824896e-11   
7.980035e-22        1000000

Some super-simple C++ code to implement this is at:
  https://www.av8n.com/physics/osc.c