BBC BASIC for Windows
Programming >> Graphics and Games >> Verlet integration
http://bb4w.conforums.com/index.cgi?board=graphics&action=display&num=1450546617
Verlet integration
Post by David Williams on Dec 19th, 2015, 4:36pm

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 )