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