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

#### STxAxTIC ##### 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.         iterations = 0
51.         PAINT (1, 1), _RGBA(0, 0, 0, 60)
52.         GOSUB drawaxes
53.
54.     ' Phase portrait
55.     CALL cpset((q1temp * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (p1temp * (2 * scalebig) + 100), Yellow)
56.     CALL cpset((q1 * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (p1 * (2 * scalebig) + 100), White)
57.     CALL cpset((q2temp * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (p2temp * (2 * scalebig) + 100), Red)
58.     CALL cpset((q2 * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (p2 * (2 * scalebig) + 100), White)
59.
60.     ' Position portrait
61.     CALL cpset((q1temp * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (q2temp * (2 * scalebig) - 100), Blue)
62.     CALL cpset((q1 * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (q2 * (2 * scalebig) - 100), White)
63.
64.     ' Pendulum
65.     CALL cline(200, 0, 400, 0, White)
66.     CALL cline(300, 0, 300 + 100 * SIN(q1temp), -100 * COS(q1temp), Black)
67.     CALL cline(300, 0, 300 + 100 * SIN(q1), -100 * COS(q1), Blue)
68.
69.     iterations = iterations + 1
70.
71.
72.
73. drawaxes:
74. ' Axis for q1 plot
75. COLOR White
76. LOCATE 2, 3: PRINT "Generalized"
77. LOCATE 3, 3: PRINT "Coordinates"
78. CALL cline(-_WIDTH / 2 + 140, 160, -_WIDTH / 2 + 400, 160, Gray)
79. COLOR White: LOCATE 5, 3: PRINT "Position q1(t)"
80. COLOR Gray: LOCATE 6, 3: PRINT "q1(0) ="; q10
81. ' Axis for p1 plot
82. CALL cline(-_WIDTH / 2 + 140, 60, -_WIDTH / 2 + 400, 60, Gray)
83. COLOR White: LOCATE 11, 3: PRINT "Momentum p1(t)"
84. COLOR Gray: LOCATE 12, 3: PRINT "p1(0) ="; p10
85. ' Axis for q2 plot
86. CALL cline(-_WIDTH / 2 + 140, -60, -_WIDTH / 2 + 400, -60, Gray)
87. COLOR White: LOCATE 18, 3: PRINT "Position q2(t)"
88. COLOR Gray: LOCATE 19, 3: PRINT "q2(0) ="; q20
89. ' Axis for p2 plot
90. CALL cline(-_WIDTH / 2 + 140, -160, -_WIDTH / 2 + 400, -160, Gray)
91. COLOR White: LOCATE 25, 3: PRINT "Momentum p2(t)"
92. COLOR Gray: LOCATE 26, 3: PRINT "p2(0) ="; p20
93. ' Axes for q & p plots
94. COLOR White
95. LOCATE 2, 60: PRINT "Phase and"
96. LOCATE 3, 58: PRINT "Position Plots"
97. LOCATE 4, 66: PRINT "p"
98. LOCATE 10, 76: PRINT "q"
99. CALL cline((190 + (320 - _WIDTH / 2)), 80 + 100, (190 + (320 - _WIDTH / 2)), -80 + 100, Gray)
100. CALL cline((190 + (320 - _WIDTH / 2)) - 100, 100, (190 + (320 - _WIDTH / 2)) + 100, 100, Gray)
101. LOCATE 17, 66: PRINT "q2"
102. LOCATE 23, 75: PRINT "q1"
103. CALL cline((190 + (320 - _WIDTH / 2)), -80 - 100, (190 + (320 - _WIDTH / 2)), 80 - 100, Gray)
104. CALL cline((190 + (320 - _WIDTH / 2)) - 100, -100, (190 + (320 - _WIDTH / 2)) + 100, -100, Gray)
105.
106. SUB cline (x1, y1, x2, y2, col AS _UNSIGNED LONG)
107.     LINE (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2)-(_WIDTH / 2 + x2, -y2 + _HEIGHT / 2), col
108.
109. SUB cpset (x1, y1, col AS _UNSIGNED LONG)
110.     PSET (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2), col
111.
« Last Edit: February 13, 2020, 06:42:19 PM by STxAxTIC »
An ounce of theory outweighs a pound of code.

#### OldMoses

• Seasoned Forum Regular
• Posts: 261 ##### 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

#### Qwerkey ##### 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

#### STxAxTIC ##### 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.

#### Qwerkey ##### 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.

#### STxAxTIC ##### Re: Pendulum Motion - Can you fix this code?
« Reply #5 on: February 14, 2020, 07:29:22 AM »

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.

#### Pete ##### 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

#### _vince ##### 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.         _PUTIMAGE, p1
47.
48.         LINE (sw/2, sh/2)-STEP(r*COS(a1-pi/2), -r*SIN(a1-pi/2))
49.         CIRCLE STEP(0,0),5
50.         LINE -STEP(r*COS(a2-pi/2), -r*SIN(a2-pi/2))
51.         CIRCLE STEP(0,0),5
52.
53.         _LIMIT 300
54.

#### STxAxTIC ##### 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.

#### Qwerkey ##### 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???

#### STxAxTIC ##### 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.

#### _vince ##### 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.
31.         LINE -(sw/2 + 100*x, sh/2 - 100*xd)
32.
33.