I've just been following a short series of YouTube-based tutorials on 'Verlet integration' from a channel called 'Coding Math'. Very good.
I adapted the code I was copying from the tutorial to attempt in an attempt to create a double pendulum simulation just so I could watch its chaotic behaviour! (Tony Tooth's program which I've downloaded more than once over the years uses a different method.)
The problem is that the physics system in the simulation listed below seems to leak a lot of energy rather quickly, and I'm not sure why. I know that a similar thing occurs with the Box2D physics engine (but the energy leak seems a lot slower!). It may just be down to a shortcoming of Verlet integration (and perhaps Euler and Runge-Kutta as well?), but if anyone can suggest how to prevent (or drastly reduce) the energy leakage, then please let me know!
EDIT: Comment from someone who tried to simulate a double pendulum using Box2D: "This is a failed experiment at a double pendulum. Box2D cannot simulate an ideal pendulum. Numerical integration causes excessive damping, which quickly makes the pendulum uninteresting."
Code: REM 'Double Pendulum' simulation using Verlet integration
REM Note: energy 'leaks' from the system quite rapidly
REM Based on the "Verlet Integration Tutorial" by Coding Math:
REM http://www.youtube.com/watch?v=3HjO_RGIjCU
REM Change to *FLOAT 64 if not using BB4W v6
*FLOAT 80
MODE 8 : OFF
ON ERROR OSCLI "REFRESH ON" : CLS : REPORT : PRINT ERL : END
WinW% = @vdu%!208
WinH% = @vdu%!212
nPts% = 3
nSticks% = 2
DIM points{(nPts%-1) pinned%, x#, y#, oldx#, oldy#, mass# }
DIM sticks{(nSticks%-1) p0%, p1%, hidden%, length# }
gravity# = -0.5
REM fixed pivot point
points{(0)}.x# = 320
points{(0)}.y# = 500
points{(0)}.oldx# = points{(0)}.x#
points{(0)}.oldy# = points{(0)}.y#
points{(0)}.mass# = 1.0
points{(0)}.pinned% = TRUE
REM pendulum bob 1
points{(1)}.x# = 320
points{(1)}.y# = 230
points{(1)}.oldx# = points{(1)}.x#
points{(1)}.oldy# = points{(1)}.y#
points{(1)}.mass# = 1.0
points{(1)}.pinned% = FALSE
REM pendulum bob 2
points{(2)}.x# = 200
points{(2)}.y# = 420
points{(2)}.oldx# = points{(2)}.x#
points{(2)}.oldy# = points{(2)}.y#
points{(2)}.mass# = 1.0
points{(2)}.pinned% = FALSE
sticks{(0)}.p0% = 0
sticks{(0)}.p1% = 1
sticks{(0)}.length# = FNstickLength#( 0 )
sticks{(0)}.hidden% = FALSE
sticks{(1)}.p0% = 1
sticks{(1)}.p1% = 2
sticks{(1)}.length# = FNstickLength#( 1 )
sticks{(1)}.hidden% = FALSE
*REFRESH OFF
REPEAT
CLS
PROCupdate
*REFRESH
WAIT 1
UNTIL FALSE
END
DEF PROCupdate
PROCupdatePoints
FOR I% = 1 TO 5
PROCupdateSticks
NEXT I%
PROCrenderPoints
PROCrenderSticks
ENDPROC
DEF PROCupdatePoints
LOCAL I%, vx#, vy#
FOR I% = 0 TO nPts%-1
IF NOT points{(I%)}.pinned% THEN
vx# = points{(I%)}.x# - points{(I%)}.oldx#
vy# = points{(I%)}.y# - points{(I%)}.oldy#
points{(I%)}.oldx# = points{(I%)}.x#
points{(I%)}.oldy# = points{(I%)}.y#
points{(I%)}.x# += vx#
points{(I%)}.y# += vy#
points{(I%)}.y# += gravity# * points{(I%)}.mass#
ENDIF
NEXT I%
ENDPROC
DEF PROCupdateSticks
LOCAL I%, dx#, dy#, distance#, difference#, percent#, offsetX#, offsetY#
FOR I% = 0 TO nSticks%-1
dx# = 1.0 * (points{(sticks{(I%)}.p1%)}.x# - points{(sticks{(I%)}.p0%)}.x#)
dy# = 1.0 * (points{(sticks{(I%)}.p1%)}.y# - points{(sticks{(I%)}.p0%)}.y#)
distance# = 1.0 * SQR(dx#^2 + dy#^2)
difference# = sticks{(I%)}.length# - distance#
percent# = difference# / distance# / 2.0
offsetX# = dx# * percent#
offsetY# = dy# * percent#
IF NOT points{( sticks{(I%)}.p0% )}.pinned% THEN
points{( sticks{(I%)}.p0% )}.x# -= offsetX#
points{( sticks{(I%)}.p0% )}.y# -= offsetY#
ENDIF
IF NOT points{( sticks{(I%)}.p1% )}.pinned% THEN
points{( sticks{(I%)}.p1% )}.x# += offsetX#
points{( sticks{(I%)}.p1% )}.y# += offsetY#
ENDIF
NEXT I%
ENDPROC
DEF PROCrenderPoints
LOCAL I%
FOR I% = 0 TO nPts%-1
CIRCLE 2*points{(I%)}.x#, 2*points{(I%)}.y#, 2*4
NEXT I%
ENDPROC
DEF PROCrenderSticks
LOCAL I%, x0%, y0%, x1%, y1%
FOR I% = 0 TO nSticks%-1
IF NOT sticks{(I%)}.hidden% THEN
x0% = points{( sticks{(I%)}.p0% )}.x#
y0% = points{( sticks{(I%)}.p0% )}.y#
x1% = points{( sticks{(I%)}.p1% )}.x#
y1% = points{( sticks{(I%)}.p1% )}.y#
LINE 2*x0%, 2*y0%, 2*x1%, 2*y1%
ENDIF
NEXT I%
ENDPROC
DEF FNstickLength#( S% )
= 1.0 * SQR( (points{( sticks{(S%)}.p0% )}.x# - points{( sticks{(S%)}.p1% )}.x#)^2 + \
\ (points{( sticks{(S%)}.p0% )}.y# - points{( sticks{(S%)}.p1% )}.y#)^2 )
Logged