### Author Topic: Physics Simulations  (Read 119 times)

#### STxAxTIC

• Library Staff
• Forum Resident
• Posts: 692
• Savage.
##### Physics Simulations
« on: February 14, 2020, 09:19:42 PM »
Hello all,

Most know I've been busy hammering down lots of stuff about computation. Pretty soon that'll be done, but to keep your whistles wet I'm posting a beefier splay of demos that examine the do's and dont's when it comes to setting up iterative equations. The program itself is dead-simple really. Most of you have already done something like what you see here.

Particle in r^-1 central potential
This applies to two-body problems in gravity or other 1/r^2 forces like electricity. There are a few boring modes, with the exciting case being (c). In the bottom-right graph, the blue ellipse is what a real planet would do, and the simulation accurately puts the sun at a focus of the ellipse (not the center).

Particle in r^-2 central potential
What if the force of gravity was different, like 1/r^3 instead of 1/r^2? Well, you'd never have been born, because orbits don't exist, so you'd have no planet. Here we show a planet trying to orbit, but it flies off into space anyway. It remains circular a while, but there's no sweet spot you can find - even on paper. (I mean analytically.)

Mass on a spring
This is a dead-simple pair of demos that show the pitfalls of Euler's method. You can crank dt way down to hide the effect.

Two equal masses connected by three springs
This one is cool. Run mode (c).

Pendulum
The remaining modes are all about the pendulum. Option 8 refers to what is for called the "Modified Euler" method (which will be, for that problem at least, a synonym for "symplectic" later.) This is what Qwerkey did. The near-critical cases are cool to watch.

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.
12. COLOR White
13. PRINT " Type a problem number and press ENTER."
14. PRINT " PROBLEM"
15. PRINT " 1) Particle in r^-1 central potential"
16. PRINT "    a) Perturbed motion ......................... Stable, Symplectic"
17. PRINT "    b) Elliptical motion ........................ Stable, Symplectic"
18. PRINT "    c) Eccentric elliptical motion .............. Stable, Symplectic"
19. PRINT " 2) Particle in r^-2 central potential"
20. PRINT "    a) Attempted circular motion ................ Unstable, Symplectic"
21. PRINT " 3) Mass on a spring"
22. PRINT "    a) Initial Momentum ..........*Incorrect*.... Unstable, Euler"
23. PRINT "    b) Initial Displacement ......*Incorrect*.... Unstable, Euler"
24. PRINT " 4) Two equal masses connected by three springs"
25. PRINT "    a) Symmetric mode ........................... Stable, Symplectic"
26. PRINT "    b) Antisymmetric mode ....................... Stable, Symplectic"
27. PRINT "    c) Perturbed motion ......................... Stable, Symplectic"
28. PRINT " 5) Plane pendulum (Part I)"
29. PRINT "    a) Initial Momentum ..........*Incorrect*.... Unstable, Euler"
30. PRINT " 6) Plane pendulum (Part III)"
31. PRINT "    a) Sub-critical case ........................ Stable, RK4"
32. PRINT "    b) Over-critical case ....................... Stable, RK4"
33. PRINT " 7) Plane pendulum (Part III)"
34. PRINT "    a) Initial Displacement ..................... Stable, Symplectic"
35. PRINT "    b) Sub-critical case ........................ Stable, Symplectic"
36. PRINT "    c) Over-critical case ....................... Stable, Symplectic"
37. PRINT " 8) Plane pendulum (Part IV)"
38. PRINT "    a) Initial Displacement ..................... Stable, Mixed Euler"; ""
39. INPUT " Enter a problem (ex. 1a): ", problem\$
40. problem\$ = LTRIM\$(RTRIM\$(LCASE\$(problem\$)))
41.
42. ' Initial conditions
43. q10 = 0
44. p10 = 0
45. q20 = 0
46. p20 = 0
47. SELECT CASE VAL(LEFT\$(problem\$, 1))
48.     CASE 1
49.         dt = .003
50.         cyclic = 1
51.         delayfactor = 1000
52.         scalebig = 20
53.         scalesmall = dt * 5
54.         m1 = 1
55.         m2 = 1
56.         g = 1
57.         SELECT CASE (RIGHT\$(problem\$, 1))
58.             CASE "a"
59.                 q10 = 1: p10 = -.2: q20 = 0: p20 = 1
60.             CASE "b"
61.                 q10 = .5: p10 = 1: q20 = -.5: p20 = 1
62.             CASE "c"
63.                 q10 = 1: p10 = .5: q20 = 0: p20 = 1.15
64.     CASE 2
65.         dt = .001
66.         cyclic = 0
67.         delayfactor = 1000
68.         scalebig = 20
69.         scalesmall = dt * 10
70.         m1 = 1
71.         m2 = 1
72.         g = 1
73.         SELECT CASE (RIGHT\$(problem\$, 1))
74.             CASE "a"
75.                 q10 = 1: p10 = 0: q20 = 0: p20 = 1.22437
76.     CASE 3
77.         dt = .03
78.         cyclic = 0
79.         delayfactor = 1000
80.         scalebig = 10
81.         scalesmall = dt * 7.5
82.         m = 1
83.         k = 1
84.         SELECT CASE (RIGHT\$(problem\$, 1))
85.             CASE "a"
86.                 q10 = 0: p10 = 2
87.             CASE "b"
88.                 q10 = 2: p10 = 0
89.     CASE 4
90.         dt = .003
91.         cyclic = 1
92.         delayfactor = 50000
93.         scalebig = 12
94.         scalesmall = dt * 7.5
95.         m = 1
96.         k = 1
97.         SELECT CASE (RIGHT\$(problem\$, 1))
98.             CASE "a"
99.                 q10 = 0: p10 = 2: q20 = 0: p20 = 2
100.             CASE "b"
101.                 q10 = 2: p10 = 0: q20 = -2: p20 = 0
102.             CASE "c"
103.                 q10 = 2: p10 = 0: q20 = 0: p20 = 0
104.     CASE 5
105.         dt = .03
106.         cyclic = 1
107.         delayfactor = 500000
108.         scalebig = 10
109.         scalesmall = dt * 5
110.         m = 1
111.         g = 1
112.         l = 1
113.         SELECT CASE (RIGHT\$(problem\$, 1))
114.             CASE "a"
115.                 q10 = 0: p10 = 1.5: scalebig = 10
116.     CASE 6
117.         dt = .03
118.         cyclic = 1
119.         delayfactor = 500000
120.         scalebig = 10
121.         scalesmall = dt * 5
122.         m = 1
123.         g = 1
124.         l = 1
125.         SELECT CASE (RIGHT\$(problem\$, 1))
126.             CASE "a"
127.                 q10 = 0: p10 = 2.19
128.             CASE "b"
129.                 q10 = 0: p10 = 2.25
130.     CASE 7, 8
131.         dt = .03
132.         cyclic = 1
133.         delayfactor = 500000
134.         scalebig = 10
135.         scalesmall = dt * 5
136.         m = 1
137.         g = 1
138.         l = 1
139.         SELECT CASE (RIGHT\$(problem\$, 1))
140.             CASE "a"
141.                 q10 = 0: p10 = -1.5
142.             CASE "b"
143.                 q10 = 0: p10 = 1.999
144.             CASE "c"
145.                 q10 = 0: p10 = 2.001
147.
148. GOSUB drawaxes
149.
150. ' Main loop.
151. iterations = 0
152. q1 = q10
153. p1 = p10
154. q2 = q20
155. p2 = p20
156.
157.
158.     FOR thedelay = 0 TO delayfactor: NEXT
159.
160.     q1temp = q1
161.     p1temp = p1
162.     q2temp = q2
163.     p2temp = p2
164.
165.     SELECT CASE VAL(LEFT\$(problem\$, 1))
166.         CASE 1
167.             ' Particle in r^-1 central potential - Symplectic integrator
168.             q1 = q1temp + dt * (p1temp / m2)
169.             q2 = q2temp + dt * (p2temp / m2)
170.             p1 = p1temp - dt * g * m1 * m2 * (q1 / ((q1 ^ 2 + q2 ^ 2) ^ (3 / 2)))
171.             p2 = p2temp - dt * g * m1 * m2 * (q2 / ((q1 ^ 2 + q2 ^ 2) ^ (3 / 2)))
172.         CASE 2
173.             ' Particle in r^-2 central potential - Symplectic integrator
174.             q1 = q1temp + dt * (p1temp / m2)
175.             q2 = q2temp + dt * (p2temp / m2)
176.             p1 = p1temp - dt * g * m1 * m2 * ((3 / 2) * q1 / ((q1 ^ 2 + q2 ^ 2) ^ (5 / 2)))
177.             p2 = p2temp - dt * g * m1 * m2 * ((3 / 2) * q2 / ((q1 ^ 2 + q2 ^ 2) ^ (5 / 2)))
178.         CASE 3
179.             ' Mass on a spring - Forward Euler method
180.             q1 = q1temp + (p1temp / m) * dt
181.             p1 = p1temp - (q1temp * k) * dt
182.             ' Mass on a spring - Backward Euler method
183.             q2 = (q2temp + (p2temp / m) * dt) / (1 + dt ^ 2)
184.             p2 = (p2temp - (q2temp * k) * dt) / (1 + dt ^ 2)
185.         CASE 4
186.             ' Two equal masses connected by three springs - Symplectic integrator
187.             q1 = q1temp + m * (p1temp) * dt
188.             p1 = p1temp - dt * k * (2 * (q1temp + m * (p1temp) * dt) - (q2temp + m * (p2temp) * dt))
189.             q2 = q2temp + m * (p2temp) * dt
190.             p2 = p2temp - dt * k * (2 * (q2temp + m * (p2temp) * dt) - (q1temp + m * (p1temp) * dt))
191.         CASE 5
192.             ' Plane pendulum - Forward Euler method
193.             q1 = q1temp + (p1temp / m) * dt
194.             p1 = p1temp - (g / l) * SIN(q1temp) * dt
195.         CASE 6
196.             ' Plane pendulum - Runge-Kutta 4th
197.             k1w = -(g / l) * SIN(q1temp)
198.             k1t = p1temp
199.             w2 = p1temp + k1w * dt / 2
200.             t2 = q1temp + k1t * dt / 2
201.             k2w = -(g / l) * SIN(t2)
202.             k2t = w2
203.             w3 = p1temp + k2w * dt / 2
204.             t3 = q1temp + k2t * dt / 2
205.             k3w = -(g / l) * SIN(t3)
206.             k3t = w3
207.             w4 = p1temp + k3w * dt
208.             t4 = q1temp + k3t * dt
209.             k4w = -(g / l) * SIN(t4)
210.             dwdt = (k1w + 2 * k2w + 2 * k3w + k4w) / 6
211.             dtdt = (k1t + 2 * k2t + 2 * k3t + k4t) / 6
212.             p1 = p1temp + dwdt * dt
213.             q1 = q1temp + dtdt * dt
214.         CASE 7
215.             ' Plane pendulum - Symplectic integrator
216.             q1 = q1temp + dt * (p1temp / m)
217.             p1 = p1temp - dt * (g / l) * (SIN(q1))
218.         CASE 8
219.             ' Plane pendulum - Modified Euler method
220.             q1 = q1temp + (p1temp / m) * dt
221.             p1 = p1temp - (g / l) * SIN(q1) * dt
222.
223.     x = (iterations * scalesmall - 180)
224.     IF (x <= 80) THEN
225.         x = x + (320 - _WIDTH / 2)
226.         CALL cpset(x, (q1 * scalebig + 160), Yellow) ' q1 plot
227.         CALL cpset(x, (p1 * scalebig + 60), Yellow) ' p1 plot
228.         CALL cpset(x, (q2 * scalebig - 60), Red) ' q2 plot
229.         CALL cpset(x, (p2 * scalebig - 160), Red) ' p2 plot
230.         IF (cyclic = 1) THEN
231.             iterations = 0
232.             PAINT (1, 1), _RGBA(0, 0, 0, 100)
233.             GOSUB drawaxes
234.
235.     ' Phase portrait
236.     CALL cpset((q1temp * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (p1temp * (2 * scalebig) + 100), Yellow)
237.     CALL cpset((q1 * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (p1 * (2 * scalebig) + 100), Gray)
238.     CALL cpset((q2temp * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (p2temp * (2 * scalebig) + 100), Red)
239.     CALL cpset((q2 * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (p2 * (2 * scalebig) + 100), White)
240.
241.     ' Position portrait
242.     CALL cpset((q1temp * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (q2temp * (2 * scalebig) - 100), Blue)
243.     CALL cpset((q1 * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (q2 * (2 * scalebig) - 100), White)
244.
245.     ' System portrait - Requires 2x wide window.
246.     SELECT CASE VAL(LEFT\$(problem\$, 1))
247.         CASE 4
248.             CALL cline(160, 0, 220 + 20 * q1, 0, Green)
249.             CALL cline(220 + 20 * q1, 0, 380 + 20 * q2, 0, Yellow)
250.             CALL cline(380 + 20 * q2, 0, 440, 0, Green)
251.             CALL ccircle(220 + 20 * q1temp, 0, 15, Black)
252.             CALL ccircle(220 + 20 * q1, 0, 15, Blue)
253.             CALL ccircle(380 + 20 * q2temp, 0, 15, Black)
254.             CALL ccircle(380 + 20 * q2, 0, 15, Red)
255.         CASE 5, 6, 7, 8
256.             CALL cline(200, 0, 400, 0, White)
257.             CALL cline(300, 0, 300 + 100 * SIN(q1temp), -100 * COS(q1temp), Black)
258.             CALL cline(300, 0, 300 + 100 * SIN(q1), -100 * COS(q1), Blue)
259.
260.     iterations = iterations + 1
261.
263.
264.
265. drawaxes:
266. ' Axis for q1 plot
267. COLOR White
268. LOCATE 2, 3: PRINT "Generalized"
269. LOCATE 3, 3: PRINT "Coordinates"
270. CALL cline(-_WIDTH / 2 + 140, 160, -_WIDTH / 2 + 400, 160, Gray)
271. COLOR White: LOCATE 5, 3: PRINT "Position q1(t)"
272. COLOR Gray: LOCATE 6, 3: PRINT "q1(0) ="; q10
273. ' Axis for p1 plot
274. CALL cline(-_WIDTH / 2 + 140, 60, -_WIDTH / 2 + 400, 60, Gray)
275. COLOR White: LOCATE 11, 3: PRINT "Momentum p1(t)"
276. COLOR Gray: LOCATE 12, 3: PRINT "p1(0) ="; p10
277. ' Axis for q2 plot
278. CALL cline(-_WIDTH / 2 + 140, -60, -_WIDTH / 2 + 400, -60, Gray)
279. COLOR White: LOCATE 18, 3: PRINT "Position q2(t)"
280. COLOR Gray: LOCATE 19, 3: PRINT "q2(0) ="; q20
281. ' Axis for p2 plot
282. CALL cline(-_WIDTH / 2 + 140, -160, -_WIDTH / 2 + 400, -160, Gray)
283. COLOR White: LOCATE 25, 3: PRINT "Momentum p2(t)"
284. COLOR Gray: LOCATE 26, 3: PRINT "p2(0) ="; p20
285. ' Axes for q & p plots
286. COLOR White
287. LOCATE 2, 60: PRINT "Phase and"
288. LOCATE 3, 58: PRINT "Position Plots"
289. LOCATE 4, 66: PRINT "p"
290. LOCATE 10, 76: PRINT "q"
291. CALL cline((190 + (320 - _WIDTH / 2)), 80 + 100, (190 + (320 - _WIDTH / 2)), -80 + 100, Gray)
292. CALL cline((190 + (320 - _WIDTH / 2)) - 100, 100, (190 + (320 - _WIDTH / 2)) + 100, 100, Gray)
293. LOCATE 17, 66: PRINT "q2"
294. LOCATE 23, 75: PRINT "q1"
295. CALL cline((190 + (320 - _WIDTH / 2)), -80 - 100, (190 + (320 - _WIDTH / 2)), 80 - 100, Gray)
296. CALL cline((190 + (320 - _WIDTH / 2)) - 100, -100, (190 + (320 - _WIDTH / 2)) + 100, -100, Gray)
297.
298. SUB cline (x1, y1, x2, y2, col AS _UNSIGNED LONG)
299.     LINE (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2)-(_WIDTH / 2 + x2, -y2 + _HEIGHT / 2), col
300.
301. SUB ccircle (x1, y1, rad, col AS _UNSIGNED LONG)
302.     CIRCLE (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2), rad, col
303.
304. SUB cpset (x1, y1, col AS _UNSIGNED LONG)
305.     PSET (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2), col
« Last Edit: February 14, 2020, 09:23:15 PM by STxAxTIC »
An ounce of theory outweighs a pound of code.