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