Author Topic: Physics Simulations  (Read 119 times)

Offline STxAxTIC

  • Library Staff
  • Forum Resident
  • Posts: 692
  • Savage.
    • Domum
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.  
  11. showmenu:
  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.         END SELECT
  65.     CASE 2
  66.         dt = .001
  67.         cyclic = 0
  68.         delayfactor = 1000
  69.         scalebig = 20
  70.         scalesmall = dt * 10
  71.         m1 = 1
  72.         m2 = 1
  73.         g = 1
  74.         SELECT CASE (RIGHT$(problem$, 1))
  75.             CASE "a"
  76.                 q10 = 1: p10 = 0: q20 = 0: p20 = 1.22437
  77.         END SELECT
  78.     CASE 3
  79.         dt = .03
  80.         cyclic = 0
  81.         delayfactor = 1000
  82.         scalebig = 10
  83.         scalesmall = dt * 7.5
  84.         m = 1
  85.         k = 1
  86.         SELECT CASE (RIGHT$(problem$, 1))
  87.             CASE "a"
  88.                 q10 = 0: p10 = 2
  89.             CASE "b"
  90.                 q10 = 2: p10 = 0
  91.         END SELECT
  92.     CASE 4
  93.         dt = .003
  94.         cyclic = 1
  95.         delayfactor = 50000
  96.         scalebig = 12
  97.         scalesmall = dt * 7.5
  98.         m = 1
  99.         k = 1
  100.         SELECT CASE (RIGHT$(problem$, 1))
  101.             CASE "a"
  102.                 q10 = 0: p10 = 2: q20 = 0: p20 = 2
  103.             CASE "b"
  104.                 q10 = 2: p10 = 0: q20 = -2: p20 = 0
  105.             CASE "c"
  106.                 q10 = 2: p10 = 0: q20 = 0: p20 = 0
  107.         END SELECT
  108.     CASE 5
  109.         dt = .03
  110.         cyclic = 1
  111.         delayfactor = 500000
  112.         scalebig = 10
  113.         scalesmall = dt * 5
  114.         m = 1
  115.         g = 1
  116.         l = 1
  117.         SELECT CASE (RIGHT$(problem$, 1))
  118.             CASE "a"
  119.                 q10 = 0: p10 = 1.5: scalebig = 10
  120.         END SELECT
  121.     CASE 6
  122.         dt = .03
  123.         cyclic = 1
  124.         delayfactor = 500000
  125.         scalebig = 10
  126.         scalesmall = dt * 5
  127.         m = 1
  128.         g = 1
  129.         l = 1
  130.         SELECT CASE (RIGHT$(problem$, 1))
  131.             CASE "a"
  132.                 q10 = 0: p10 = 2.19
  133.             CASE "b"
  134.                 q10 = 0: p10 = 2.25
  135.         END SELECT
  136.     CASE 7, 8
  137.         dt = .03
  138.         cyclic = 1
  139.         delayfactor = 500000
  140.         scalebig = 10
  141.         scalesmall = dt * 5
  142.         m = 1
  143.         g = 1
  144.         l = 1
  145.         SELECT CASE (RIGHT$(problem$, 1))
  146.             CASE "a"
  147.                 q10 = 0: p10 = -1.5
  148.             CASE "b"
  149.                 q10 = 0: p10 = 1.999
  150.             CASE "c"
  151.                 q10 = 0: p10 = 2.001
  152.         END SELECT
  153.     CASE ELSE
  154.         GOTO showmenu
  155.  
  156. GOSUB drawaxes
  157.  
  158. ' Main loop.
  159. iterations = 0
  160. q1 = q10
  161. p1 = p10
  162. q2 = q20
  163. p2 = p20
  164.  
  165.  
  166.     FOR thedelay = 0 TO delayfactor: NEXT
  167.  
  168.     q1temp = q1
  169.     p1temp = p1
  170.     q2temp = q2
  171.     p2temp = p2
  172.  
  173.     SELECT CASE VAL(LEFT$(problem$, 1))
  174.         CASE 1
  175.             ' Particle in r^-1 central potential - Symplectic integrator
  176.             q1 = q1temp + dt * (p1temp / m2)
  177.             q2 = q2temp + dt * (p2temp / m2)
  178.             p1 = p1temp - dt * g * m1 * m2 * (q1 / ((q1 ^ 2 + q2 ^ 2) ^ (3 / 2)))
  179.             p2 = p2temp - dt * g * m1 * m2 * (q2 / ((q1 ^ 2 + q2 ^ 2) ^ (3 / 2)))
  180.         CASE 2
  181.             ' Particle in r^-2 central potential - Symplectic integrator
  182.             q1 = q1temp + dt * (p1temp / m2)
  183.             q2 = q2temp + dt * (p2temp / m2)
  184.             p1 = p1temp - dt * g * m1 * m2 * ((3 / 2) * q1 / ((q1 ^ 2 + q2 ^ 2) ^ (5 / 2)))
  185.             p2 = p2temp - dt * g * m1 * m2 * ((3 / 2) * q2 / ((q1 ^ 2 + q2 ^ 2) ^ (5 / 2)))
  186.         CASE 3
  187.             ' Mass on a spring - Forward Euler method
  188.             q1 = q1temp + (p1temp / m) * dt
  189.             p1 = p1temp - (q1temp * k) * dt
  190.             ' Mass on a spring - Backward Euler method
  191.             q2 = (q2temp + (p2temp / m) * dt) / (1 + dt ^ 2)
  192.             p2 = (p2temp - (q2temp * k) * dt) / (1 + dt ^ 2)
  193.         CASE 4
  194.             ' Two equal masses connected by three springs - Symplectic integrator
  195.             q1 = q1temp + m * (p1temp) * dt
  196.             p1 = p1temp - dt * k * (2 * (q1temp + m * (p1temp) * dt) - (q2temp + m * (p2temp) * dt))
  197.             q2 = q2temp + m * (p2temp) * dt
  198.             p2 = p2temp - dt * k * (2 * (q2temp + m * (p2temp) * dt) - (q1temp + m * (p1temp) * dt))
  199.         CASE 5
  200.             ' Plane pendulum - Forward Euler method
  201.             q1 = q1temp + (p1temp / m) * dt
  202.             p1 = p1temp - (g / l) * SIN(q1temp) * dt
  203.         CASE 6
  204.             ' Plane pendulum - Runge-Kutta 4th
  205.             k1w = -(g / l) * SIN(q1temp)
  206.             k1t = p1temp
  207.             w2 = p1temp + k1w * dt / 2
  208.             t2 = q1temp + k1t * dt / 2
  209.             k2w = -(g / l) * SIN(t2)
  210.             k2t = w2
  211.             w3 = p1temp + k2w * dt / 2
  212.             t3 = q1temp + k2t * dt / 2
  213.             k3w = -(g / l) * SIN(t3)
  214.             k3t = w3
  215.             w4 = p1temp + k3w * dt
  216.             t4 = q1temp + k3t * dt
  217.             k4w = -(g / l) * SIN(t4)
  218.             dwdt = (k1w + 2 * k2w + 2 * k3w + k4w) / 6
  219.             dtdt = (k1t + 2 * k2t + 2 * k3t + k4t) / 6
  220.             p1 = p1temp + dwdt * dt
  221.             q1 = q1temp + dtdt * dt
  222.         CASE 7
  223.             ' Plane pendulum - Symplectic integrator
  224.             q1 = q1temp + dt * (p1temp / m)
  225.             p1 = p1temp - dt * (g / l) * (SIN(q1))
  226.         CASE 8
  227.             ' Plane pendulum - Modified Euler method
  228.             q1 = q1temp + (p1temp / m) * dt
  229.             p1 = p1temp - (g / l) * SIN(q1) * dt
  230.     END SELECT
  231.  
  232.     x = (iterations * scalesmall - 180)
  233.     IF (x <= 80) THEN
  234.         x = x + (320 - _WIDTH / 2)
  235.         CALL cpset(x, (q1 * scalebig + 160), Yellow) ' q1 plot
  236.         CALL cpset(x, (p1 * scalebig + 60), Yellow) ' p1 plot
  237.         CALL cpset(x, (q2 * scalebig - 60), Red) ' q2 plot
  238.         CALL cpset(x, (p2 * scalebig - 160), Red) ' p2 plot
  239.     ELSE
  240.         IF (cyclic = 1) THEN
  241.             iterations = 0
  242.             PAINT (1, 1), _RGBA(0, 0, 0, 100)
  243.             GOSUB drawaxes
  244.         ELSE
  245.             EXIT DO
  246.         END IF
  247.     END IF
  248.  
  249.     ' Phase portrait
  250.     CALL cpset((q1temp * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (p1temp * (2 * scalebig) + 100), Yellow)
  251.     CALL cpset((q1 * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (p1 * (2 * scalebig) + 100), Gray)
  252.     CALL cpset((q2temp * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (p2temp * (2 * scalebig) + 100), Red)
  253.     CALL cpset((q2 * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (p2 * (2 * scalebig) + 100), White)
  254.  
  255.     ' Position portrait
  256.     CALL cpset((q1temp * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (q2temp * (2 * scalebig) - 100), Blue)
  257.     CALL cpset((q1 * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (q2 * (2 * scalebig) - 100), White)
  258.  
  259.     ' System portrait - Requires 2x wide window.
  260.     SELECT CASE VAL(LEFT$(problem$, 1))
  261.         CASE 4
  262.             CALL cline(160, 0, 220 + 20 * q1, 0, Green)
  263.             CALL cline(220 + 20 * q1, 0, 380 + 20 * q2, 0, Yellow)
  264.             CALL cline(380 + 20 * q2, 0, 440, 0, Green)
  265.             CALL ccircle(220 + 20 * q1temp, 0, 15, Black)
  266.             CALL ccircle(220 + 20 * q1, 0, 15, Blue)
  267.             CALL ccircle(380 + 20 * q2temp, 0, 15, Black)
  268.             CALL ccircle(380 + 20 * q2, 0, 15, Red)
  269.         CASE 5, 6, 7, 8
  270.             CALL cline(200, 0, 400, 0, White)
  271.             CALL cline(300, 0, 300 + 100 * SIN(q1temp), -100 * COS(q1temp), Black)
  272.             CALL cline(300, 0, 300 + 100 * SIN(q1), -100 * COS(q1), Blue)
  273.     END SELECT
  274.  
  275.     iterations = iterations + 1
  276.  
  277. GOTO showmenu
  278.  
  279.  
  280. drawaxes:
  281. ' Axis for q1 plot
  282. COLOR White
  283. LOCATE 2, 3: PRINT "Generalized"
  284. LOCATE 3, 3: PRINT "Coordinates"
  285. CALL cline(-_WIDTH / 2 + 140, 160, -_WIDTH / 2 + 400, 160, Gray)
  286. COLOR White: LOCATE 5, 3: PRINT "Position q1(t)"
  287. COLOR Gray: LOCATE 6, 3: PRINT "q1(0) ="; q10
  288. ' Axis for p1 plot
  289. CALL cline(-_WIDTH / 2 + 140, 60, -_WIDTH / 2 + 400, 60, Gray)
  290. COLOR White: LOCATE 11, 3: PRINT "Momentum p1(t)"
  291. COLOR Gray: LOCATE 12, 3: PRINT "p1(0) ="; p10
  292. ' Axis for q2 plot
  293. CALL cline(-_WIDTH / 2 + 140, -60, -_WIDTH / 2 + 400, -60, Gray)
  294. COLOR White: LOCATE 18, 3: PRINT "Position q2(t)"
  295. COLOR Gray: LOCATE 19, 3: PRINT "q2(0) ="; q20
  296. ' Axis for p2 plot
  297. CALL cline(-_WIDTH / 2 + 140, -160, -_WIDTH / 2 + 400, -160, Gray)
  298. COLOR White: LOCATE 25, 3: PRINT "Momentum p2(t)"
  299. COLOR Gray: LOCATE 26, 3: PRINT "p2(0) ="; p20
  300. ' Axes for q & p plots
  301. COLOR White
  302. LOCATE 2, 60: PRINT "Phase and"
  303. LOCATE 3, 58: PRINT "Position Plots"
  304. LOCATE 4, 66: PRINT "p"
  305. LOCATE 10, 76: PRINT "q"
  306. CALL cline((190 + (320 - _WIDTH / 2)), 80 + 100, (190 + (320 - _WIDTH / 2)), -80 + 100, Gray)
  307. CALL cline((190 + (320 - _WIDTH / 2)) - 100, 100, (190 + (320 - _WIDTH / 2)) + 100, 100, Gray)
  308. LOCATE 17, 66: PRINT "q2"
  309. LOCATE 23, 75: PRINT "q1"
  310. CALL cline((190 + (320 - _WIDTH / 2)), -80 - 100, (190 + (320 - _WIDTH / 2)), 80 - 100, Gray)
  311. CALL cline((190 + (320 - _WIDTH / 2)) - 100, -100, (190 + (320 - _WIDTH / 2)) + 100, -100, Gray)
  312.  
  313. SUB cline (x1, y1, x2, y2, col AS _UNSIGNED LONG)
  314.     LINE (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2)-(_WIDTH / 2 + x2, -y2 + _HEIGHT / 2), col
  315.  
  316. SUB ccircle (x1, y1, rad, col AS _UNSIGNED LONG)
  317.     CIRCLE (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2), rad, col
  318.  
  319. SUB cpset (x1, y1, col AS _UNSIGNED LONG)
  320.     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.