Chronology Current Month Current Thread Current Date
[Year List] [Month List (current year)] [Date Index] [Thread Index] [Thread Prev] [Thread Next] [Date Prev] [Date Next]

Re: [Phys-L] Help w/ Euler Cromer algo.



On 2014, Sep 23, , at 10:59, Bill Nettles <bnettles@uu.edu> wrote:

> BC:
> I'm trying to reproduce your original problem and am finding that the 
> residual (thetamax*cos(wt)-theta_calculated(t)) depends on thetamax.   What 
> was your theta max to get such small residuals (~e-9)?  If I use a thetamax 
> of 0.1 radians, I get a maximum residual of about 3.1*e-4 over 5 periods and 
> it doesn't grow quickly. That's with 1 ms time step.  When I cut the time to 
> 0.1 ms, the residual drops to 3.1*e-5.  That's a first order effect which is 
> what I would expect.  My residual is not growing quickly. It grows from 3 to 
> 3.7 (*e-4) over 50 periods. At 100 periods (100,000 time steps) it's barely 
> over 4.
> 
> I'm using iPython notebook to compute this.  Looks pretty stable to me.


Wow!  thanks for doing this.

Do you have Kaleidagraph?

I use TrueBASIC on an old G4 macbook. and then Kal.  

I’m about to leave town, so will reply;  return.  I’ll, perforce, put pics of 
output, etc. on my site [ http://www.cleyet.org ]  link to  in reply.

However, here’s a boiler plate of what I’m currently using:   ! = interpreter 
(or compiler) ignore. PRINT  “ … “  = print as string in output.

In most runs I use 0.1 ms (deltat).   To delineraize the oscillator:    theta=> 
sine(theta)  [in Q= line — coulda done DO line]

The g line(s) was for when I was investigating the “lunar effect”   Remember 
the long conversations on adiabatic invariance?


=====================
PRINT  "boiler plate for lin. Damped Simple Pendulum w/ Sync. type impulse and 
g change"
PRINT "this version test of lin nothing else"
PRINT  "initial displacement angle in radians 5 deg. ~=88mrad"
OPTION ANGLE radians
LET theta = 88e-3
!INPUT theta
LET thetadot = 0
LET t = 0
PRINT "time increment (delta t)"
!let deltat = 1e-3
INPUT deltat
!print deltat
LET g = 9.82
LET l = 0.2487435058              ! w / above ~ a one second period pendulum 
PRINT "damping coefficient if not inputted is 0.0126 (Q=500)"
PRINT "is a one s pendulum (accuracy ~ 10 figs.)"
!LET r = 0.0126
INPUT r
PRINT "drive amplitude"
INPUT d
PRINT "g increment (using? 1e-5 if not inputted g inc. = 0)"
INPUT deltag
!LET deltag = 0
PRINT "number of steps (1k steps / sec.) if deltat is 1ms"
INPUT stop
PRINT "numbermod =10 w/ deltat 1ms => 100 prints/s; 100=>10p/s"
INPUT numbermod
PRINT "  t(seconds)   Theta(radians)"  !   speed      g"
PRINT t,theta                     !,thetadot,g
SUB p
    PRINT t,theta                 !,thetadot,g
END SUB
SUB step
    !IF t >=5 then LET dc = r ELSE let dc = 0  !sets when damped
    IF theta < 0.05 AND theta > 0.0 AND thetadot > 0.0 THEN LET da = d ELSE LET 
da = 0.0      !sets when driven
    LET Q = - (g/l) * theta  - r*thetadot +da    ! Q = d^2 (theta) / dt^2
    LET thetadot = Q * deltat + thetadot
    LET theta = theta + deltat * thetadot

    LET t = t + deltat            ! increments time
END SUB
FOR N = 1 to stop
    CALL step
    IF mod (N, numbermod) = 0 and t> 0.0 then CALL p  !mod80:1e 400:0.2e
    !IF abs(thetadot) <= 0.2 then CALL p             ! find a max.
    !IF (t => 2) and t< 48 then LET g = g + deltag    ! delta g added 
continuously beginning at three seconds to 17.
NEXT N
END
=================

bc