Forced Damped Pendulum

Here are some questions for you to think about and things to do.

  1. What are the real equations of motion?  What are the units on these various constants?
  2. Can you reduce the number of independent parameters, by making things dimensionless?
  
  3. Browse through the mathematica help files for the functions below, i.e.  Do, ListPlot, Interpolation, NDSolve,...
  
  4. What happens qualitatively to the solution when the friction K gets smaller?
       Or when the forcing term B gets bigger?  
       Or when the forcing frequence ω is near the natural frequency?  (What is the natural frequencey?)
  
  5. For some combination of parameters, the motion settles down into a steady pattern.
      For others it doesn't.  What is the dividing line between the two?  Can you make some kind of 
      plot (or color-coded plot, even if only by hand) of what the various regions of "parameter phase space" are?
      
  6. I said in class that phase space plots don't cross lines.  But the one below does - how come?
  
  7. Can you alter the program (or write one of your own) that finds the points where the forcing term
      has the same phase each time?  What does that look in (x,v) phase space?
      
  8.  This system exhibits both regular and chaotic motion.  Yeah, I know we haven't defined those things
        rigorously yet, but can you find examples of parameters where each is going on?  
        
  9.  Do the initial starting conditions matter for the long-term behavior?  What do two trajectories that
       start near each other do over time?  (This is part of the answer to 8.)
       
 10. How does the accuracy of solving this vary with the numerical method?  
        I'm using either a built-in one (which you hope the Mathematica folks did right), or the simplest Euler method.
        Others to try would be predictor-corrector, Runge-Kutta, ...
	
11. Can you get  numbers from some C code into Mathematica ?  (See the ReadList function, 2.11.7 in the Mathematica Book)

Equation of motion :    (((∂^( 2)) x(t))/(∂ (t^2)))= -A sin(x)  -  K (( ∂ x)/∂) t + 
             B sin(ω t)

accel[x_,  t_,  v_]  :=  
      - A Sin[x] + B Sin[ω t] -  K v;

Definition of the parameters

A  =  1;  B  = 0.9;   K =  0.3;  ω  =  1.8; 
    tMaxi  = 50;

Euler's method of numerically solving it : constant time step, save all the results in Mathematica lists.  (Ugly syntax for doing so, I'm afraid.)

\(x\  = \ 0.0; v\  = \ 1.0; 
dt\  = \ 0.02; 
xpoints\  = \ {}; 
    vpoints = {}; phipoints = {}; 

    Do[\ 
		\ v\  = \ v\  + \ accel[x, t, v]\ dt; 
		\ 
      x\  = \ x\  + \ v\ dt; 
		xpoints = Append[xpoints, {t, x}]; 
		
      \(vpoints = Append[vpoints, {t, v}];\)
	, \ {t, 0, tMaxi, dt}
	]
      \)ListPlot[xpoints,  PlotJoined -> True]

And here's the phase space plot of the last calculation.

\(xp\  = \ Interpolation[xpoints]; 
vp\  = \ Interpolation[vpoints]; 

    ParametricPlot[{xp[t], vp[t]}, {t, 0, tMaxi}]\)

Or, use the built-in NDSolve : Numerical DiffEQ Solver.  (If it doesn't hang on you...)

\(anser\  = \ 
      NDSolve[\ {\(\(y'\)'\)[t] == accel[y[t], t, \(y'\)[t]], y[0] == 0, \ 
          \(y'\)[0] == 1}, y, {t, 0, tMaxi}]; 

    Plot[Evaluate[y[t]\  /. \ anser], {t, 0, tMaxi}]; \)