PHYSICS 115/242

    Example solution in C, for Qu. 6, HW 3

This question requires a bit more thought. We first of all need to
decide when the particle has completed one oscillation.  I did this by
determining the time when the momentum changes sign. This is half the
period so I multiplied this time by 2. (I did a linear extrapolation
between the last two times to determine more precisely the time at which
the momentum vanished).

The other issue is to estimate whether the time step is small enough
that reasonable accuracy has been achieved. I did this by checking
whether the value of x after a half period is equal to -x_0 (where x_0
is the initial displacement, i.e. the amplitude of the oscillation) to
within some tolerance, which I took to be 10^(-3).  The source code
below incorporates these features and uses 2nd order Runge Kutta (RK2)
to integrate the motion.

Note:
1. I checked the code for the simple harmonic oscillator, by setting
f(x) = -x and checking that the period was (close to) 2 \pi for any
initial displacement.

2. During exectution, the program asks me for the values of the step
size and the starting value for x. In this way I can change these
parameters without having to edit the code and recompile. Notice the
structure of the scanf command.

===========================================================================
Source Code
===========================================================================

#include "math.h"
#include "stdio.h"

double f(double x)                      // define the force
{
    return -x*x*x ;
}

main()
{
    double x, f(), p, dt, t, h, k1x, k2x, k1p, k2p, x0, pold, thalf, eps;
    int half_orbit;

    printf (" x0 = ? h = ? " );
    scanf("%lf %lf", &x0, &h);
    printf (" x0 = %10.5f  h = %10.5f \n", x0, h) ;

    eps = 0.001;
    half_orbit = 0;
    x = x0;
    p = 0;
    t = 0;

    while (half_orbit == 0 )
    {
        pold = p;
        k1x =   p;               // 2nd order Runge-Kutta
        k1p = f(x);
        k2x =   p + h * k1p;
        k2p = f(x + h * k1x);
        x = x + 0.5 * h * (k1x + k2x);
        p = p + 0.5 * h * (k1p + k2p);
        t = t + h;
        if (p > 0) 	// We've done half an orbit
        {
            half_orbit = 1;
            dt = h * p / ( p - pold);
            thalf = t - dt;                //linear interpolation to determine
                                           // time where p went through zero
            if (fabs(x+x0) > eps) 
            {
                printf ( "error in x too great; need to reduce h \n");
            }
        }
    }
    printf (" period = %10.5f \n", 2 * thalf );
}

===========================================================================
Output: Code
===========================================================================
 ~/courses/115/C => a.out
 x0 = ? h = ? 0.1 0.1
 x0 =    0.10000  h =    0.10000 
 period =   74.16066 
 ~/courses/115/C => a.out
 x0 = ? h = ? 1 0.05
 x0 =    1.00000  h =    0.05000 
 period =    7.41035 
 ~/courses/115/C => a.out
 x0 = ? h = ? 10 0.005
 x0 =   10.00000  h =    0.00500 
error in x too great; need to reduce h 
 period =    0.74104 
 ~/courses/115/C => a.out
 x0 = ? h = ? 10 0.002
 x0 =   10.00000  h =    0.00200 
 period =    0.74154 

These results give the output for x0 = 0.1, 1, and 10. Note that a very
small stepsize is needed to get the desired accuray for x0 = 10.

===========================================================================
Comments:
===========================================================================

For the simple harmonic oscillator the force and (hence the velocity)
increase as the amplitude increases. Hence, although the particle has
further to go, it goes more quickly and these two effects compensate in
the determination of the period. For the x^4 potential the force
increases VERY fast as the amplitude increases, so the increase in
velocity is greater than the increase in distance, so the period goes
down (like 1/amplitude) as the amplitude increases. This can be shown by
dimensional analysis.