Author Topic: Pendulum Motion - Can you fix this code?  (Read 201 times)

Offline STxAxTIC

  • Library Staff
  • Forum Addict
  • Posts: 560
  • Savage.
    • Domum
Pendulum Motion - Can you fix this code?
« on: February 13, 2020, 06:32:43 PM »
Hello all,

I've taken several pains to prepare a range of QB64 demos, all in the scheme of fleshing out methods for calculating things exact(enough)ly. This goes from numbers, trajectories, differential equations, and so on, leaving almost no stone unturned. Pretty soon I'll be finished with a hefty PDF on all this, in where I will formally argue against what I can imagine some of you are thinking right now: "I don't need no stinkin' methods, I just tweak my program til it looks right". I won't truly argue with these people (this time). Your methods work for you, mine work for me, live and let live. If you don't mind though, I'm going to continue the imaginary argument we're having for the sake of motivating the issues. Just pretend.

All that said, this post is NOT even about the final product. For now, I want to make people familiar with the graphs and notation I'm using. When you run the code below, and see the screenshots, you will see a handful of graphs. Now don't just tuck your head into your turtleneck - these are easy. On the right, you see a pendulum swinging in gravity - that's the *real* system we're simulating. On the left I track (as q1) the angle from vertical displacement of the pendulum, and also I track the angular velocity (as p1). There are other graphs that will be super obvious when you watch. Just like regular displacement and regular velocity, q1 and p1 determine the motion. Together these are, as you may have guessed, sinusoidal (swingy). Look at the sine-looking graphs on the left. All easy.

Except this code has a problem. (On purpose.) See how the pendulum keeps on swinging higher without being asked? That's an error in the code. I gotta be honest, pretty much all code by non-science students, and 90% of people who "should" know better, code the same error, only they unknowingly hide it using very small time steps, so the error looks small... for a while... This demo uses a big time step on purpose to exaggerate the effect, which is why I had to slow it down like crazy. Tweak that as needed.

The part that updates the motion is

Code: QB64: [Select]
  1.     q1 = q1temp + (p1temp / m) * dt
  2.     p1 = p1temp - (g / l) * SIN(q1temp) * dt

... which is as simple as it gets for a pendulum. The combination p1temp/m is playing the role of velocity, and (g/l)*sin(q1temp) is proportional to the downward force. Just like a ball flying through the air, just has a SIN() in there.

So you're thinking "Okay, I'll just always use the equations I already like, but will keep a small time step and never suffer this effect".... *BZZZT* Nope, such code will still be wrong eventually, as you're completely unable to simulate periodic motion or any long-lasting trajectory. Plus, you want your program to be fast, right? If you have the math right, you *can* take big time steps and have the motion still be right. Otherwise you're stuck going slow if you want the motion to be right.

Pretty soon I'll post a satisfactory answer to all this, but for now: Who can make the pendulum swing right? If you get it by blind luck, great. If you get it by deduction, even better. If you come up with a method that works for all problems, you certainly don't need me. The solution has to still look correct in gravity. No windshield wipers.

Code: QB64: [Select]
  1. CONST Black = _RGB32(0, 0, 0)
  2. CONST Blue = _RGB32(0, 0, 255)
  3. CONST Gray = _RGB32(128, 128, 128)
  4. CONST Green = _RGB32(0, 128, 0)
  5. CONST Red = _RGB32(255, 0, 0)
  6. CONST White = _RGB32(255, 255, 255)
  7. CONST Yellow = _RGB32(255, 255, 0)
  8.  
  9. SCREEN _NEWIMAGE(1280, 480, 32)
  10.  
  11. ' Initial conditions
  12. dt = .03
  13. delayfactor = 500000
  14.  
  15. m = 1 ' Mass
  16. g = 1 ' Gravity constant
  17. l = 1 ' Pendulum length
  18. q10 = 0 ' Initial angle (vertical)
  19. p10 = 1.5 ' Initial momentum (a kick to the right)
  20. scalebig = 10: scalesmall = 5 * dt
  21. GOSUB drawaxes
  22.  
  23. ' Main loop.
  24. iterations = 0
  25. q1 = q10
  26. p1 = p10
  27. q2 = q20
  28. p2 = p20
  29.  
  30.     FOR thedelay = 0 TO delayfactor: NEXT
  31.  
  32.     q1temp = q1
  33.     p1temp = p1
  34.     q2temp = q2
  35.     p2temp = p2
  36.  
  37.     ' Plane pendulum calculations
  38.     ''' What's wrong with these?
  39.     q1 = q1temp + (p1temp / m) * dt
  40.     p1 = p1temp - (g / l) * SIN(q1temp) * dt
  41.     ''' These aren't quite right...
  42.  
  43.     x = (iterations * scalesmall - 180)
  44.     IF (x <= 80) THEN
  45.         x = x + (320 - _WIDTH / 2)
  46.         CALL cpset(x, (q1 * scalebig + 160), Yellow) ' q1 plot
  47.         CALL cpset(x, (p1 * scalebig + 60), Yellow) ' p1 plot
  48.         CALL cpset(x, (q2 * scalebig - 60), Red) ' q2 plot
  49.         CALL cpset(x, (p2 * scalebig - 160), Red) ' p2 plot
  50.     ELSE
  51.         iterations = 0
  52.         PAINT (1, 1), _RGBA(0, 0, 0, 60)
  53.         GOSUB drawaxes
  54.     END IF
  55.  
  56.     ' Phase portrait
  57.     CALL cpset((q1temp * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (p1temp * (2 * scalebig) + 100), Yellow)
  58.     CALL cpset((q1 * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (p1 * (2 * scalebig) + 100), White)
  59.     CALL cpset((q2temp * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (p2temp * (2 * scalebig) + 100), Red)
  60.     CALL cpset((q2 * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (p2 * (2 * scalebig) + 100), White)
  61.  
  62.     ' Position portrait
  63.     CALL cpset((q1temp * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (q2temp * (2 * scalebig) - 100), Blue)
  64.     CALL cpset((q1 * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (q2 * (2 * scalebig) - 100), White)
  65.  
  66.     ' Pendulum
  67.     CALL cline(200, 0, 400, 0, White)
  68.     CALL cline(300, 0, 300 + 100 * SIN(q1temp), -100 * COS(q1temp), Black)
  69.     CALL cline(300, 0, 300 + 100 * SIN(q1), -100 * COS(q1), Blue)
  70.  
  71.     iterations = iterations + 1
  72.  
  73.  
  74.  
  75. drawaxes:
  76. ' Axis for q1 plot
  77. COLOR White
  78. LOCATE 2, 3: PRINT "Generalized"
  79. LOCATE 3, 3: PRINT "Coordinates"
  80. CALL cline(-_WIDTH / 2 + 140, 160, -_WIDTH / 2 + 400, 160, Gray)
  81. COLOR White: LOCATE 5, 3: PRINT "Position q1(t)"
  82. COLOR Gray: LOCATE 6, 3: PRINT "q1(0) ="; q10
  83. ' Axis for p1 plot
  84. CALL cline(-_WIDTH / 2 + 140, 60, -_WIDTH / 2 + 400, 60, Gray)
  85. COLOR White: LOCATE 11, 3: PRINT "Momentum p1(t)"
  86. COLOR Gray: LOCATE 12, 3: PRINT "p1(0) ="; p10
  87. ' Axis for q2 plot
  88. CALL cline(-_WIDTH / 2 + 140, -60, -_WIDTH / 2 + 400, -60, Gray)
  89. COLOR White: LOCATE 18, 3: PRINT "Position q2(t)"
  90. COLOR Gray: LOCATE 19, 3: PRINT "q2(0) ="; q20
  91. ' Axis for p2 plot
  92. CALL cline(-_WIDTH / 2 + 140, -160, -_WIDTH / 2 + 400, -160, Gray)
  93. COLOR White: LOCATE 25, 3: PRINT "Momentum p2(t)"
  94. COLOR Gray: LOCATE 26, 3: PRINT "p2(0) ="; p20
  95. ' Axes for q & p plots
  96. COLOR White
  97. LOCATE 2, 60: PRINT "Phase and"
  98. LOCATE 3, 58: PRINT "Position Plots"
  99. LOCATE 4, 66: PRINT "p"
  100. LOCATE 10, 76: PRINT "q"
  101. CALL cline((190 + (320 - _WIDTH / 2)), 80 + 100, (190 + (320 - _WIDTH / 2)), -80 + 100, Gray)
  102. CALL cline((190 + (320 - _WIDTH / 2)) - 100, 100, (190 + (320 - _WIDTH / 2)) + 100, 100, Gray)
  103. LOCATE 17, 66: PRINT "q2"
  104. LOCATE 23, 75: PRINT "q1"
  105. CALL cline((190 + (320 - _WIDTH / 2)), -80 - 100, (190 + (320 - _WIDTH / 2)), 80 - 100, Gray)
  106. CALL cline((190 + (320 - _WIDTH / 2)) - 100, -100, (190 + (320 - _WIDTH / 2)) + 100, -100, Gray)
  107.  
  108. SUB cline (x1, y1, x2, y2, col AS _UNSIGNED LONG)
  109.     LINE (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2)-(_WIDTH / 2 + x2, -y2 + _HEIGHT / 2), col
  110.  
  111. SUB cpset (x1, y1, col AS _UNSIGNED LONG)
  112.     PSET (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2), col
  113.  
« Last Edit: February 13, 2020, 06:42:19 PM by STxAxTIC »
An ounce of theory outweighs a pound of code.

Offline OldMoses

  • Forum Regular
  • Posts: 240
Re: Pendulum Motion - Can you fix this code?
« Reply #1 on: February 13, 2020, 06:44:09 PM »
I'm just blown away by how much action you can put in less than 90 lines of working code...
Andy

Offline Qwerkey

  • Forum Addict
  • Posts: 507
Re: Pendulum Motion - Can you fix this code?
« Reply #2 on: February 14, 2020, 05:28:15 AM »
STxAxTIC, my Cuckoo Clock program used pendulum (gravity controlled) maths with no problems (as far as I was aware) of the position drifting off.  The simulation maths is not that difficult, I thought.

https://www.qb64.org/forum/index.php?topic=1115.msg103115#msg103115


Offline STxAxTIC

  • Library Staff
  • Forum Addict
  • Posts: 560
  • Savage.
    • Domum
Re: Pendulum Motion - Can you fix this code?
« Reply #3 on: February 14, 2020, 06:58:21 AM »
Hello Qwerkey,

I just looked at your code, and yes - that motion in fact IS correct! Cool. So this means you fall into one of three categories:

Quote
Who can make the pendulum swing right?

(A) If you get it by blind luck, great.
(B) If you get it by deduction, even better.
(C) If you come up with a method that works for all problems, you certainly don't need me.

Without spoiling the answer, you know which of group (A), (B), or (C) you belong to. Please tell us if you like!
An ounce of theory outweighs a pound of code.

Offline Qwerkey

  • Forum Addict
  • Posts: 507
Re: Pendulum Motion - Can you fix this code?
« Reply #4 on: February 14, 2020, 07:19:35 AM »
STxAxTIC, I think that you may have missed a fourth category in which I place myself(!):

(D) A supremely-gifted coder who has never needed the least amount of help from Steve McNeill and whose grasp of the most difficult mathematics is legendary and who knows no bounds of self-aggrandisement.


Offline STxAxTIC

  • Library Staff
  • Forum Addict
  • Posts: 560
  • Savage.
    • Domum
Re: Pendulum Motion - Can you fix this code?
« Reply #5 on: February 14, 2020, 07:29:22 AM »
Lol, username checks out!

Actually if one stares are the two essential lines of your code versus the two essential lines above, the difference stands out. It's literally a change of ______. Too bad it only works for the pendulum though, and maybe other O(2) problems at best.
An ounce of theory outweighs a pound of code.

Offline Pete

  • Forum Addict
  • Posts: 1695
  • Cuz I sez so, varmint!
Re: Pendulum Motion - Can you fix this code?
« Reply #6 on: February 14, 2020, 11:37:37 AM »
If there's going to be a D added, there needs to be an E as well, as in...

E) NONE OF THE ABOVE!!!

How's it swinging, Bill?

Pete :D

Offline _vince

  • Forum Regular
  • Posts: 167
Re: Pendulum Motion - Can you fix this code?
« Reply #7 on: February 15, 2020, 03:46:04 AM »
i'll leave this

Code: QB64: [Select]
  1. DEFSNG a-z
  2. pi = 3.141593
  3. r = 100
  4. a1 = -pi/2
  5. a2 = -pi/2
  6. h=0.001
  7. m=1000
  8. g=1000
  9. sw = 640
  10. sh = 480
  11.  
  12. p1 = _NEWIMAGE(sw,sh,12)
  13. SCREEN _NEWIMAGE(sw,sh,12)
  14.  
  15.         a1p = (6/(m*r^2))*(2*m*v1-3*COS(a1-a2)*m*v2)/(16-9*(COS(a1-a2))^2)
  16.         a1k1 = a1p
  17.         a1k2 = a1p + h*a1k1/2
  18.         a1k3 = a1p + h*a1k2/2
  19.         a1k4 = a1p + h*a1k3
  20.         a2p = (6/(m*r^2))*(8*m*v2-3*COS(a1-a2)*m*v1)/(16-9*(COS(a1-a2))^2)
  21.         a2k1 = a2p
  22.         a2k2 = a2p + h*a2k1/2
  23.         a2k3 = a2p + h*a2k2/2
  24.         a2k4 = a2p + h*a2k3
  25.         na1=a1+h*(a1k1+2*a1k2+2*a1k3+a1k4)/6
  26.         na2=a2+h*(a2k1+2*a2k2+2*a2k3+a2k4)/6
  27.         v1p = -0.5*r^2*(a1p*a2p*SIN(a1-a2)+3*g*SIN(a1)/r)
  28.         v1k1 = v1p
  29.         v1k2 = v1p + h*v1k1/2
  30.         v1k3 = v1p + h*v1k2/2
  31.         v1k4 = v1p + h*v1k3
  32.         v2p = -0.5*r^2*(-a1p*a2p*SIN(a1-a2)+g*SIN(a2)/r)
  33.         v2k1 = v2p
  34.         v2k2 = v2p + h*v2k1/2
  35.         v2k3 = v2p + h*v2k2/2
  36.         v2k4 = v2p + h*v2k3
  37.         nv1=v1+h*(v1k1+2*v1k2+2*v1k3+v1k4)/6
  38.         nv2=v2+h*(v2k1+2*v2k2+2*v2k3+v2k4)/6
  39.         a1=na1
  40.         a2=na2
  41.         v1=nv1
  42.         v2=nv2
  43.  
  44.         _DEST p1
  45.         PSET (sw/2+r*COS(a1-pi/2)+r*COS(a2-pi/2), sh/2-r*SIN(a1-pi/2)-r*SIN(a2-pi/2)),12
  46.         _DEST 0
  47.         _PUTIMAGE, p1
  48.  
  49.         LINE (sw/2, sh/2)-STEP(r*COS(a1-pi/2), -r*SIN(a1-pi/2))
  50.         CIRCLE STEP(0,0),5
  51.         LINE -STEP(r*COS(a2-pi/2), -r*SIN(a2-pi/2))
  52.         CIRCLE STEP(0,0),5
  53.  
  54.         _LIMIT 300
  55.         _DISPLAY
  56.  

Offline STxAxTIC

  • Library Staff
  • Forum Addict
  • Posts: 560
  • Savage.
    • Domum
Re: Pendulum Motion - Can you fix this code?
« Reply #8 on: February 15, 2020, 06:24:27 AM »
Hi_vince,

Ah, the good old double pendulum solved by the fourth-order Runge-Kitta technique. *Someone* has been reading physics books!

Nice work. Surely you can fix the original broken code and answer the challenge by swapping in the entire four-order apparatus, but can do you do it in one line? hehe
An ounce of theory outweighs a pound of code.

Offline Qwerkey

  • Forum Addict
  • Posts: 507
Re: Pendulum Motion - Can you fix this code?
« Reply #9 on: February 15, 2020, 07:50:31 AM »
... by swapping in the entire four-order apparatus...

Oh, crikey!! What is he talking about???

Offline STxAxTIC

  • Library Staff
  • Forum Addict
  • Posts: 560
  • Savage.
    • Domum
Re: Pendulum Motion - Can you fix this code?
« Reply #10 on: February 15, 2020, 08:33:26 AM »
Oh, crikey!! What is he talking about???


I'm talking about mode 6 here:

https://www.qb64.org/forum/index.php?topic=2186.0

See this form repeated in _vince's code?

Code: QB64: [Select]
  1.             ' Plane pendulum - Runge-Kutta 4th
  2.             k1w = -(g / l) * SIN(q1temp)
  3.             k1t = p1temp
  4.             w2 = p1temp + k1w * dt / 2
  5.             t2 = q1temp + k1t * dt / 2
  6.             k2w = -(g / l) * SIN(t2)
  7.             k2t = w2
  8.             w3 = p1temp + k2w * dt / 2
  9.             t3 = q1temp + k2t * dt / 2
  10.             k3w = -(g / l) * SIN(t3)
  11.             k3t = w3
  12.             w4 = p1temp + k3w * dt
  13.             t4 = q1temp + k3t * dt
  14.             k4w = -(g / l) * SIN(t4)
  15.             dwdt = (k1w + 2 * k2w + 2 * k3w + k4w) / 6
  16.             dtdt = (k1t + 2 * k2t + 2 * k3t + k4t) / 6
  17.             p1 = p1temp + dwdt * dt
  18.             q1 = q1temp + dtdt * dt
An ounce of theory outweighs a pound of code.

Offline _vince

  • Forum Regular
  • Posts: 167
Re: Pendulum Motion - Can you fix this code?
« Reply #11 on: February 16, 2020, 08:22:38 PM »
The double pendulum is a well known 'chaotic' system. Here's an excerpt you'd appreciate

Code: QB64: [Select]
  1. sw = 800
  2. sh = 600
  3.  
  4. 'screenres sw, sh, 32
  5. 'windowtitle "kicked harmonic oscillator phase-space trajectory"
  6. SCREEN _NEWIMAGE(sw,sh,32)
  7.  
  8.  
  9. a = 1
  10. b = 1
  11. w = 1
  12. tt = 5
  13. k = 0
  14.  
  15. t = 0
  16. x = a*COS(w*t) + b*SIN(w*t)
  17. xd = -w*a*SIN(w*t) + w*b*COS(w*t)
  18. PSET (sw/2 + 100*x, sh/2 - 100*xd)
  19.  
  20. FOR t = 0 TO 100 STEP 0.1
  21.         x = a*COS(w*t) + b*SIN(w*t)
  22.         xd = -w*a*SIN(w*t) + w*b*COS(w*t)
  23.  
  24.         IF INT(t/tt) > k THEN
  25.                 a = a - SIN(w*t)/w
  26.                 b = b + COS(w*t)/w
  27.  
  28.                 k = INT(t/tt)
  29.                 'pset (sw/2 + 100*x, sh/2 - 100*xd)
  30.         END IF
  31.  
  32.         LINE -(sw/2 + 100*x, sh/2 - 100*xd)
  33.  
  34.