hp41programs

Alea Double Pendulum for the HP-41

Overview

O
| \
| \
| \ L1
| \
| \
| \M
| /
| / L2
| /
| N
|
|
|
x

-With

µ₁ = ( Ox , OM ) & µ₂ = ( Ox , MN )
L₁ & L₂ = lengths of rod (constant)
m₁ & m₂ = mass of pendulum
g = gravitational constant

-The differential equations are:

µ₁" = [ −g (2 m₁ + m₂) sin µ₁ − m₂ g sin(µ₁ − 2 µ₂) − 2 sin(µ₁ − µ₂) m₂ (µ'₂² L₂ + µ'₁² L₁ cos(µ₁ − µ₂)) ] / [ L₁ (2 m₁ + m₂ − m₂ cos(2 µ₁ − 2 µ₂)) ]
µ₂" = [ 2 sin(µ₁−µ₂) (µ'₁² L₁ (m₁ + m₂) + g(m₁ + m₂) cos µ₁ + µ'₂² L₂ m₂ cos(µ₁ − µ₂)) ] / [ L₂ (2 m₁ + m₂ − m₂ cos(2 µ₁ − 2 µ₂)) ]

-"DPN" employs a Runge-Kutta method of order 4:

Program Listing

Data Registers: R00 & R09-R10 & R15 to R26: temp ( • Registers are to be initialized before executing "DPN" )

• R01 = t₀ • R04 = m₁ • R06 = L₁ • R08 = g • R11 = µ₁ ( rd ) • R13 = µ'₁ ( rd/s )
• R02 = h = stepsize • R05 = m₂ • R07 = L₂ • R12 = µ₂ ( rd ) • R14 = µ'₂ ( rd/s )
• R03 = N = number of steps

Flags: /
Subroutines: /

-Line 228 is a three-byte GTO 01

01 LBL "DPN"
02 RAD
03 RCL 02
04 2
05 /
06 STO 09
07 RCL 03
08 STO 16
09 GTO 01
10 LBL 00
11 RCL 13
12 X^2
13 RCL 06
14 *
15 RCL 11
16 RCL 12
17 -
18 1
19 P-R
20 STO 10
21 X<>Y
22 STO 00
23 RDN
24 *
25 RCL 14
26 X^2
27 RCL 07
28 *
29 +
30 RCL 05
31 *
32 RCL 00
33 *
34 ST+ X

35 RCL 11
36 RCL 12
37 ST+ X
38 -
39 SIN
40 RCL 05
41 *
42 RCL 04
43 ST+ X
44 LASTX
45 +
46 RCL 11
47 SIN
48 *
49 +
50 RCL 08
51 *
52 +
53 RCL 06
54 /
55 CHS
56 STO 15
57 RCL 05
58 RCL 10
59 *
60 RCL 07
61 *
62 RCL 14
63 X^2
64 *
65 RCL 11
66 COS
67 RCL 08
68 *

69 RCL 13
70 X^2
71 RCL 06
72 *
73 +
74 RCL 04
75 RCL 05
76 +
77 *
78 +
79 RCL 00
80 ST+ X
81 *
82 RCL 07
83 /
84 RCL 04
85 ST+ X
86 RCL 05
87 +
88 RCL 11
89 RCL 12
90 -
91 ST+ X
92 COS
93 RCL 05
94 *
95 -
96 RCL 02
97 /
98 ST/ 15
99 /
100 RTN
101 LBL 01
102 RCL 11

103 STO 17
104 RCL 12
105 STO 18
106 RCL 13
107 STO 19
108 RCL 14
109 STO 20
110 XEQ 00
111 STO 22
112 2
113 /
114 ST+ 14
115 RCL 15
116 STO 21
117 LASTX
118 /
119 ST+ 13
120 LASTX
121 ST/ Z
122 /
123 RCL 19
124 +
125 RCL 09
126 *
127 ST+ 11
128 CLX
129 RCL 20
130 +
131 RCL 09
132 *
133 ST+ 12
134 XEQ 00
135 STO 24
136 2

137 /
138 RCL 20
139 +
140 STO 14
141 RCL 15
142 STO 23
143 2
144 /
145 RCL 19
146 +
147 STO 13
148 XEQ 00
149 STO 26
150 RCL 20
151 +
152 STO 14
153 RCL 15
154 STO 25
155 RCL 19
156 +
157 STO 13
158 RCL 26
159 2
160 ST/ 15
161 /
162 RCL 20
163 +
164 RCL 02
165 *
166 RCL 18
167 +
168 STO 12
169 RCL 15

170 RCL 19
171 +
172 RCL 02
173 *
174 RCL 17
175 +
176 STO 11
177 XEQ 00
178 STO 10
179 RCL 21
180 RCL 23
181 RCL 25
182 +
183 STO 00
184 +
185 ST+ 00
186 6
187 /
188 RCL 19
189 STO 13
190 +
191 RCL 02
192 *
193 RCL 17
194 +
195 STO 11
196 RCL 00
197 RCL 15
198 +
199 6
200 /
201 ST+ 13
202 RCL 22

203 RCL 24
204 RCL 26
205 +
206 STO 00
207 +
208 ST+ 00
209 6
210 /
211 RCL 20
212 STO 14
213 +
214 RCL 02
215 *
216 RCL 18
217 +
218 STO 12
219 RCL 00
220 RCL 10
221 +
222 6
223 /
224 ST+ 14
225 RCL 02
226 ST+ 01
227 DSE 16
228 GTO 01
229 RCL 01
230 SIGN
231 RCL 14
232 RCL 13
233 RCL 12
234 RCL 11
235 END

( 305 bytes / SIZE 027 )

STACK	INPUTS	OUTPUTS
T	/	µ'₂
Z	/	µ'₁
Y	/	µ₂
X	_/	µ₁
L	/	t₀ + N.h

Example: t₀ = 0 h = 0.01 N = 10

m₁ = 1 L₁ = 3 µ₁ = µ'₁ = µ'₂ = 0 µ₂ = 1 radian g = 9.81
m₂ = 2 L₂ = 4

0 STO 01 0.01 STO 02 10 STO 03 0 STO 11 STO 13 STO 14 1 STO 04 STO 12 2 STO 05 3 STO 06 4 STO 07 9.81 STO 08

XEQ "DPN" >>>> µ₁= 0.006224775 rd ---Execution time = 4mn02s---
RDN µ₂ = 0.987152548 rd
RDN µ'₁ = 0.125934133 rd/s
RDN µ'₂ = -0.257668985 rd/s

and LASTX = t = 0.1 s

-Then, R/S >>>> µ₁= 0.025775346 rd
   RDN µ₂ = 0.948159032 rd
   RDN µ'₁ = 0.269616097 rd/s
   RDN µ'₂ = -0.524685660 rd/s

and LASTX = t = 0.2 s

Reference:

[1] https://web.mit.edu/jorloff/www/chaosTalk/double-pendulum/double-pendulum-en.html