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
q1 = q1temp + (p1temp / m) * dt
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.
' Initial conditions
dt = .03
delayfactor = 500000
m = 1 ' Mass
g = 1 ' Gravity constant
l = 1 ' Pendulum length
q10 = 0 ' Initial angle (vertical)
p10 = 1.5 ' Initial momentum (a kick to the right)
scalebig = 10: scalesmall = 5 * dt
' Main loop.
iterations = 0
q1 = q10
p1 = p10
q2 = q20
p2 = p20
q1temp = q1
p1temp = p1
q2temp = q2
p2temp = p2
' Plane pendulum calculations
''' What's wrong with these?
q1 = q1temp + (p1temp / m) * dt
p1
= p1temp
- (g
/ l
) * SIN(q1temp
) * dt
''' These aren't quite right...
x = (iterations * scalesmall - 180)
CALL cpset
(x
, (q1
* scalebig
+ 160), Yellow
) ' q1 plot CALL cpset
(x
, (p1
* scalebig
+ 60), Yellow
) ' p1 plot CALL cpset
(x
, (q2
* scalebig
- 60), Red
) ' q2 plot CALL cpset
(x
, (p2
* scalebig
- 160), Red
) ' p2 plot iterations = 0
' Phase portrait
CALL cpset
((q1temp
* (2 * scalebig
) + (190 + (320 - _WIDTH / 2))), (p1temp
* (2 * scalebig
) + 100), Yellow
) CALL cpset
((q1
* (2 * scalebig
) + (190 + (320 - _WIDTH / 2))), (p1
* (2 * scalebig
) + 100), White
) CALL cpset
((q2temp
* (2 * scalebig
) + (190 + (320 - _WIDTH / 2))), (p2temp
* (2 * scalebig
) + 100), Red
) CALL cpset
((q2
* (2 * scalebig
) + (190 + (320 - _WIDTH / 2))), (p2
* (2 * scalebig
) + 100), White
)
' Position portrait
CALL cpset
((q1temp
* (2 * scalebig
) + (190 + (320 - _WIDTH / 2))), (q2temp
* (2 * scalebig
) - 100), Blue
) CALL cpset
((q1
* (2 * scalebig
) + (190 + (320 - _WIDTH / 2))), (q2
* (2 * scalebig
) - 100), White
)
' Pendulum
CALL cline
(200, 0, 400, 0, White
) CALL cline
(300, 0, 300 + 100 * SIN(q1temp
), -100 * COS(q1temp
), Black
) CALL cline
(300, 0, 300 + 100 * SIN(q1
), -100 * COS(q1
), Blue
)
iterations = iterations + 1
drawaxes:
' Axis for q1 plot
' Axis for p1 plot
' Axis for q2 plot
' Axis for p2 plot
' Axes for q & p plots
CALL cline
((190 + (320 - _WIDTH / 2)), 80 + 100, (190 + (320 - _WIDTH / 2)), -80 + 100, Gray
) CALL cline
((190 + (320 - _WIDTH / 2)) - 100, 100, (190 + (320 - _WIDTH / 2)) + 100, 100, Gray
) CALL cline
((190 + (320 - _WIDTH / 2)), -80 - 100, (190 + (320 - _WIDTH / 2)), 80 - 100, Gray
) CALL cline
((190 + (320 - _WIDTH / 2)) - 100, -100, (190 + (320 - _WIDTH / 2)) + 100, -100, Gray
)
