Text Only | Text with Attachments

QB64.org Forum

Active Forums => Programs => Topic started by: STxAxTIC on February 13, 2020, 06:32:43 pm

Title: Pendulum Motion - Can you fix this code?
Post by: STxAxTIC 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)
  11.  
  12. ' Initial conditions
  13. dt = .03
  14. delayfactor = 500000
  15.  
  16. m = 1 ' Mass
  17. g = 1 ' Gravity constant
  18. l = 1 ' Pendulum length
  19. q10 = 0 ' Initial angle (vertical)
  20. p10 = 1.5 ' Initial momentum (a kick to the right)
  21. scalebig = 10: scalesmall = 5 * dt
  22. GOSUB drawaxes
  23.  
  24. ' Main loop.
  25. iterations = 0
  26. q1 = q10
  27. p1 = p10
  28. q2 = q20
  29. p2 = p20
  30.  
  32.     FOR thedelay = 0 TO delayfactor: NEXT
  33.  
  34.     q1temp = q1
  35.     p1temp = p1
  36.     q2temp = q2
  37.     p2temp = p2
  38.  
  39.     ' Plane pendulum calculations
  40.     ''' What's wrong with these?
  41.     q1 = q1temp + (p1temp / m) * dt
  42.     p1 = p1temp - (g / l) * SIN(q1temp) * dt
  43.     ''' These aren't quite right...
  44.  
  45.     x = (iterations * scalesmall - 180)
  46.     IF (x <= 80) THEN
  47.         x = x + (320 - _WIDTH / 2)
  48.         CALL cpset(x, (q1 * scalebig + 160), Yellow) ' q1 plot
  49.         CALL cpset(x, (p1 * scalebig + 60), Yellow) ' p1 plot
  50.         CALL cpset(x, (q2 * scalebig - 60), Red) ' q2 plot
  51.         CALL cpset(x, (p2 * scalebig - 160), Red) ' p2 plot
  52.     ELSE
  53.         iterations = 0
  54.         PAINT (1, 1), _RGBA(0, 0, 0, 60)
  55.         GOSUB drawaxes
  56.     END IF
  57.  
  58.     ' Phase portrait
  59.     CALL cpset((q1temp * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (p1temp * (2 * scalebig) + 100), Yellow)
  60.     CALL cpset((q1 * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (p1 * (2 * scalebig) + 100), White)
  61.     CALL cpset((q2temp * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (p2temp * (2 * scalebig) + 100), Red)
  62.     CALL cpset((q2 * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (p2 * (2 * scalebig) + 100), White)
  63.  
  64.     ' Position portrait
  65.     CALL cpset((q1temp * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (q2temp * (2 * scalebig) - 100), Blue)
  66.     CALL cpset((q1 * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (q2 * (2 * scalebig) - 100), White)
  67.  
  68.     ' Pendulum
  69.     CALL cline(200, 0, 400, 0, White)
  70.     CALL cline(300, 0, 300 + 100 * SIN(q1temp), -100 * COS(q1temp), Black)
  71.     CALL cline(300, 0, 300 + 100 * SIN(q1), -100 * COS(q1), Blue)
  72.  
  73.     iterations = iterations + 1
  74.  
  76.  
  78.  
  79. drawaxes:
  80. ' Axis for q1 plot
  81. COLOR White
  82. LOCATE 2, 3: PRINT "Generalized"
  83. LOCATE 3, 3: PRINT "Coordinates"
  84. CALL cline(-_WIDTH / 2 + 140, 160, -_WIDTH / 2 + 400, 160, Gray)
  85. COLOR White: LOCATE 5, 3: PRINT "Position q1(t)"
  86. COLOR Gray: LOCATE 6, 3: PRINT "q1(0) ="; q10
  87. ' Axis for p1 plot
  88. CALL cline(-_WIDTH / 2 + 140, 60, -_WIDTH / 2 + 400, 60, Gray)
  89. COLOR White: LOCATE 11, 3: PRINT "Momentum p1(t)"
  90. COLOR Gray: LOCATE 12, 3: PRINT "p1(0) ="; p10
  91. ' Axis for q2 plot
  92. CALL cline(-_WIDTH / 2 + 140, -60, -_WIDTH / 2 + 400, -60, Gray)
  93. COLOR White: LOCATE 18, 3: PRINT "Position q2(t)"
  94. COLOR Gray: LOCATE 19, 3: PRINT "q2(0) ="; q20
  95. ' Axis for p2 plot
  96. CALL cline(-_WIDTH / 2 + 140, -160, -_WIDTH / 2 + 400, -160, Gray)
  97. COLOR White: LOCATE 25, 3: PRINT "Momentum p2(t)"
  98. COLOR Gray: LOCATE 26, 3: PRINT "p2(0) ="; p20
  99. ' Axes for q & p plots
  100. COLOR White
  101. LOCATE 2, 60: PRINT "Phase and"
  102. LOCATE 3, 58: PRINT "Position Plots"
  103. LOCATE 4, 66: PRINT "p"
  104. LOCATE 10, 76: PRINT "q"
  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. LOCATE 17, 66: PRINT "q2"
  108. LOCATE 23, 75: PRINT "q1"
  109. CALL cline((190 + (320 - _WIDTH / 2)), -80 - 100, (190 + (320 - _WIDTH / 2)), 80 - 100, Gray)
  110. CALL cline((190 + (320 - _WIDTH / 2)) - 100, -100, (190 + (320 - _WIDTH / 2)) + 100, -100, Gray)
  112.  
  113. SUB cline (x1, y1, x2, y2, col AS _UNSIGNED LONG)
  114.     LINE (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2)-(_WIDTH / 2 + x2, -y2 + _HEIGHT / 2), col
  116.  
  117. SUB cpset (x1, y1, col AS _UNSIGNED LONG)
  118.     PSET (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2), col
  120.  
Title: Re: Pendulum Motion - Can you fix this code?
Post by: OldMoses 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...
Title: Re: Pendulum Motion - Can you fix this code?
Post by: Qwerkey 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 (https://www.qb64.org/forum/index.php?topic=1115.msg103115#msg103115)

Title: Re: Pendulum Motion - Can you fix this code?
Post by: STxAxTIC 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!
Title: Re: Pendulum Motion - Can you fix this code?
Post by: Qwerkey 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.

Title: Re: Pendulum Motion - Can you fix this code?
Post by: STxAxTIC 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.
Title: Re: Pendulum Motion - Can you fix this code?
Post by: Pete 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
Title: Re: Pendulum Motion - Can you fix this code?
Post by: _vince 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.  
  13. p1 = _NEWIMAGE(sw,sh,12)
  14. SCREEN _NEWIMAGE(sw,sh,12)
  15.  
  17.         a1p = (6/(m*r^2))*(2*m*v1-3*COS(a1-a2)*m*v2)/(16-9*(COS(a1-a2))^2)
  18.         a1k1 = a1p
  19.         a1k2 = a1p + h*a1k1/2
  20.         a1k3 = a1p + h*a1k2/2
  21.         a1k4 = a1p + h*a1k3
  22.         a2p = (6/(m*r^2))*(8*m*v2-3*COS(a1-a2)*m*v1)/(16-9*(COS(a1-a2))^2)
  23.         a2k1 = a2p
  24.         a2k2 = a2p + h*a2k1/2
  25.         a2k3 = a2p + h*a2k2/2
  26.         a2k4 = a2p + h*a2k3
  27.         na1=a1+h*(a1k1+2*a1k2+2*a1k3+a1k4)/6
  28.         na2=a2+h*(a2k1+2*a2k2+2*a2k3+a2k4)/6
  29.         v1p = -0.5*r^2*(a1p*a2p*SIN(a1-a2)+3*g*SIN(a1)/r)
  30.         v1k1 = v1p
  31.         v1k2 = v1p + h*v1k1/2
  32.         v1k3 = v1p + h*v1k2/2
  33.         v1k4 = v1p + h*v1k3
  34.         v2p = -0.5*r^2*(-a1p*a2p*SIN(a1-a2)+g*SIN(a2)/r)
  35.         v2k1 = v2p
  36.         v2k2 = v2p + h*v2k1/2
  37.         v2k3 = v2p + h*v2k2/2
  38.         v2k4 = v2p + h*v2k3
  39.         nv1=v1+h*(v1k1+2*v1k2+2*v1k3+v1k4)/6
  40.         nv2=v2+h*(v2k1+2*v2k2+2*v2k3+v2k4)/6
  41.         a1=na1
  42.         a2=na2
  43.         v1=nv1
  44.         v2=nv2
  45.  
  46.         _DEST p1
  47.         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
  48.         _DEST 0
  49.         _PUTIMAGE, p1
  50.  
  51.         LINE (sw/2, sh/2)-STEP(r*COS(a1-pi/2), -r*SIN(a1-pi/2))
  52.         CIRCLE STEP(0,0),5
  53.         LINE -STEP(r*COS(a2-pi/2), -r*SIN(a2-pi/2))
  54.         CIRCLE STEP(0,0),5
  55.  
  56.         _LIMIT 300
  57.         _DISPLAY
  60.  
Title: Re: Pendulum Motion - Can you fix this code?
Post by: STxAxTIC 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
Title: Re: Pendulum Motion - Can you fix this code?
Post by: Qwerkey on February 15, 2020, 07:50:31 am
Quote from: STxAxTIC on February 15, 2020, 06:24:27 am
... by swapping in the entire four-order apparatus...

Oh, crikey!! What is he talking about???
Title: Re: Pendulum Motion - Can you fix this code?
Post by: STxAxTIC on February 15, 2020, 08:33:26 am
Quote from: Qwerkey on February 15, 2020, 07:50:31 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 (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
Title: Re: Pendulum Motion - Can you fix this code?
Post by: _vince 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.  
  15.  
  16. a = 1
  17. b = 1
  18. w = 1
  19. tt = 5
  20. k = 0
  21.  
  22. t = 0
  23. x = a*cos(w*t) + b*sin(w*t)
  24. xd = -w*a*sin(w*t) + w*b*cos(w*t)
  25. pset (sw/2 + 100*x, sh/2 - 100*xd)
  26.  
  27. for t = 0 to 100 step 0.1
  28.         x = a*cos(w*t) + b*sin(w*t)
  29.         xd = -w*a*sin(w*t) + w*b*cos(w*t)
  30.  
  31.         if int(t/tt) > k then
  32.                 a = a - sin(w*t)/w
  33.                 b = b + cos(w*t)/w
  34.  
  35.                 k = int(t/tt)
  36.                 'pset (sw/2 + 100*x, sh/2 - 100*xd)
  37.         end if
  38.  
  39.         line -(sw/2 + 100*x, sh/2 - 100*xd)
  41.  
  45.  
Text Only | Text with Attachments