« previous next »
Pages: [1]   Go Down

Author Topic: Physics Simulations  (Read 3186 times)

0 Members and 1 Guest are viewing this topic.

Offline STxAxTIC

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