DefDbl A-Z: DefInt I-J, N
n = 5: nn = 0
Dim x(n), xx(n), f(n), c(4, n)
nextone:
Screen 0: Cls
_FullScreen
' option here to get the x-scale NOT magnified as in popular modern monitors
GoSub graphics
GoSub startup
agn:
Cls
Line (0, -a)-(0, a), 8
ph = -a * Cos(w * x(5))
Circle (0, ph), .025, 12: Paint (0, ph), 12
X1 = Sin(xx(1))
X2 = -Cos(xx(1)) + ph
X3 = Sin(xx(1)) + Sin(xx(3))
X4 = -Cos(xx(1)) - Cos(xx(3)) + ph
Line (0, ph)-(X1, X2), 12
Circle (X1, X2), .05, 9: Paint (X1, X2), 9
Line (X1, X2)-(X3, X4), 12
Circle (X3, X4), .05, 9: Paint (X3, X4), 9
_Delay .001
GoSub Runge
t = t + h
Locate 21, 6: Print CSng(t)
a$ = InKey$: If a$ = "" Then GoTo agn Else GoTo nextone
Equations:
Q = x(3) - x(1)
c = Cos(Q)
s = Sin(Q)
P = c * s
D = 1# - m * c ^ 2#
g = (1# + a * w ^ 2# * Cos(w * x(5)))
g1 = g * Sin(x(1))
g3 = g * Sin(x(3))
x22 = x(2) ^ 2#
x42 = x(4) ^ 2#
f(1) = x(2)
f(2) = -k * x(2) + (m * (x42 * s + x22 * P + g3 * c) - g1) / D
f(3) = x(4)
f(4) = -k * x(4) + (-x22 * s - m * x42 * P - g3 + g1 * c) / D
f(5) = 1#
Locate 7, 5: Print x(1)
Locate 8, 5: Print x(2)
Locate 9, 5: Print x(3)
Locate 10, 5: Print x(4)
Locate 11, 5: Print x(5)
_Display
Return
Runge:
For i = 1 To n: x(i) = xx(i): Next
GoSub Equations
For i = 1 To n: c(1, i) = h * f(i): Next
For i = 1 To n: x(i) = xx(i) + c(1, i) / 2#: Next
GoSub Equations
For i = 1 To n: c(2, i) = h * f(i): Next
For i = 1 To n: x(i) = xx(i) + c(2, i) / 2#: Next
GoSub Equations
For i = 1 To n: c(3, i) = h * f(i): Next
For i = 1 To n: x(i) = xx(i) + c(3, i): Next
GoSub Equations
For i = 1 To n: c(4, i) = h * f(i): Next
For i = 1 To n
xx(i) = xx(i) + (c(1, i) + 2# * c(2, i) + 2# * c(3, i) + c(4, i)) / 6#
Next
Return
graphics:
Cls: Screen 9
Paint (1, 1), 9
View (220, 17)-(595, 330), 0, 14
Window (-3.6, -2.4)-(3.6, 2.8)
Locate 21, 2: Print "Time"
Return
startup:
k = .1#: m = .1#: xx(2) = 0: xx(4) = 0
'nn = 2
Select Case nn
Case O
k = .05#: m = .1#: xx(2) = 0: xx(4) = 0: xx(1) = .4: xx(3) = -.2
Print "viscous damping k=0.05"
Case 1
k = 0#: m = .1#: xx(2) = 0: xx(4) = 0: xx(1) = 0: xx(3) = 3.1428
Locate 5, 2
Print " no viscous damping "
Case 2
a = .2: w = 10: xx(1) = 3.1: xx(3) = 1 'stable yet semi-hanging
Locate 3, 1
Print "stable yet top-hanging"
Case 3
a = .1: w = 3.696#: xx(1) = .01: xx(3) = .01: m = .5# 'fast mode
Locate 3, 1
Print " fast-pump mode "
Case 4
a = .1: w = 1.55#: xx(1) = .01: xx(3) = .01: m = .5#: k = .03# 'slow mode
Locate 3, 1
Print " slow-pump mode "
Case 5
a = .1: w = 1.55#: xx(1) = .1: xx(3) = .2: m = .5# 'k=.1 too much for slow
Locate 3, 1
Print " damping k=0.1 is too "
Locate 4, 1
Print " much for slow mode "
Case 6
a = .1: w = 2: xx(1) = .32: xx(3) = .32 'goes to down-hanging
Locate 4, 1
Print " goes to down-hanging "
Case 7
a = .1: w = 2: xx(1) = 1: xx(3) = -1 'stable fast swing
Locate 3, 2
Print " stable fast swing "
Case 8
a = .25: w = 2: xx(1) = 1.5: xx(3) = -1.5 'stable fast swing
Locate 3, 2
Print " stable fast swing "
Case 9
a = .25#: w = 2#: xx(1) = 1: xx(3) = -1: x(2) = 5: x(4) = 5: k = .05 'whirling
Locate 3, 2
Print " whirling dervrish "
Case 10
m = .5#: k = .2#: xx(1) = 2.5: xx(3) = 1.9: a = .1: w = 25 'a<.3;w>1.85/a
Locate 3, 2
Print " Indian Rope Trick! "
Case 11
m = .5#: k = .2#: xx(1) = 2.6: xx(3) = -3.1: a = .1: w = 25 'a<.3;w>1.85/a
Locate 3, 3
Print " teeter-totter "
Case 12
m = .5#: k = .2#: xx(1) = 3.1: xx(3) = 3.2: a = .04: w = 25 'collapse 1
Locate 3, 4
Print " surprise! "
Case 13
m = .5#: k = .2#: xx(1) = 3.1: xx(3) = 3.2: a = .15#: w = 25 'collapse 2
Locate 3, 2
Print " this IS demo 13! "
Case 14
m = .5#: k = .2#: xx(1) = 3.05#: xx(3) = 3.25#: a = .11#: w = 25 'a<.3;w>1.85/a
Case Is > 15
Screen 0: Cls: Print "Only maths and physics types live in a LINEAR delusion"
Print " where you can uncook eggs and correct errors!": End
End Select
If nn > 1 Then
Locate 5, 4
Print " damping k=";: Print Using "#.##"; k
Else
Locate 3, 3
Print " A double pendulum"
Locate 4, 3
Print " swinging with "
End If
Locate 25, 1: Print "press any key for next demo";
'IF nn < 0 THEN LOCATE 12, 2: INPUT "a,w"; a, w
t = 0#: xx(5) = 0#
'IF nn > 0 THEN LOCATE 14, 2: INPUT "ang1,ang2"; xx(1), xx(3)
h = .005#
nn = nn + 1
Return