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.



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.
For 1 period:   0.1 ms time step -> 3.8e-8   and 1.0 ms time step -> 1.2 e-6

The differences (in my calcs) oscillate like a sine curve which says that the 
largest discrepancies are occurring near theta=0 and smallest near + and - 
thetamax.

-> -----Original Message-----
-> From: Phys-l [mailto:phys-l-bounces@www.phys-l.org] On Behalf Of
-> Bernard Cleyet
-> Sent: Tuesday, September 23, 2014 2:50 PM
-> To: Phys-L@Phys-L.org
-> Subject: 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
-> _______________________________________________
-> Forum for Physics Educators
-> Phys-l@www.phys-l.org
-> http://www.phys-l.org/mailman/listinfo/phys-l