Our goal is to simulate a mechanical system according to its Hamiltonian $\mathcal H\left(\mathbf q,\mathbf p,t\right)$. To utilize the canonical equations \begin{equation} \frac{\mathrm d\mathbf q}{\mathrm dt}= \frac{\partial\mathcal H}{\partial\mathbf p},\quad \frac{\mathrm d\mathbf p}{\mathrm dt}= -\frac{\partial\mathcal H}{\partial\mathbf p}, \label{canonical} \end{equation} we need to calculate the partial derivatives of $\mathcal H$. Here is a simple code to calculate partial derivatives. def div x0, dx, f f0 = f.(x0) n = x0.size Array.new n do |i| (f.(x0 + Vector.basis(n, i) dx) - f0) / dx end end (RGSS do not have matrix.rb, you can copy one from the attached file below.) Here x0 is a Vector, f is a call-able object as a function of vectors, dx is a small scalar which we are going to take 1e-6. Let x = Vector[q, p], and then Formula \ref{canonical} has the form \begin{equation} \frac{\mathrm d\mathbf x}{\mathrm dt}=f\left(\mathbf x\right). \end{equation} To solve this equation numerically, we need to use a famous method called the (explicit) Runge–Kutta method. def runge_kutta initial, max_t, dt, (pyramid, coefs), func (0..max_t).step(dt).reduce initial do |ret, t| $canvas.trace t, ret if $canvas coefs.zip(pyramid).each_with_object([]).sum do |(coef, row), ary| coef ary.push(func.(t, row.inner(ary) dt + ret)).last end dt + ret#p(ret) end end Note that Runge–Kutta is a family of methods. The argument (pyramid, coefs) takes one of the following, each of which is a single Runge–Kutta method. FORWARD_EULER = [[],[1]] EXPLICIT_MIDPOINT = [[],[1/2.0],[0,1]] HEUN = [[],[1],[1/2.0,1/2.0]] RALSTON = [[],[2/3.0],[1/4.0,3/4.0]] KUTTA_3RD = [[],[1/2.0],[-1,2],[1/6.0,2/3.0,1/6.0]] HEUN_3RD = [[],[1/3.0],[0,2/3.0],[1/4.0,0,3/4.0]] RALSTON_3RD = [[],[1/2.0],[0,3/4.0],[2/9.0,1/3.0,4/9.0]] SSPRK3 = [[],[1],[1/4.0,1/4.0],[1/6.0,1/6.0,2/3.0]] CLASSIC_4TH = [[],[1/2.0],[0,1/2.0],[0,0,1],[1/6.0,1/3.0,1/3.0,1/6.0]] RALSTON_4TH = [[],[0.4],[0.29697761,0.15875964],[0.21810040,-3.05096516, 3.83286476],[0.17476028, -0.55148066, 1.20553560, 0.17118478]] THREE_EIGHTH_4TH = [[],[1/3.0],[-1/3.0,1],[1,-1,1],[1/8.0,3/8.0,3/8.0,1/8.0]] Here we are going to take CLASSIC_4TH. The $canvas appearing here is an object that is going to draw the result onto the screen. Here we also need to have some patches to get it work. class Float alias ulysses20200426121236_add + def + other if zero? && [Vector, Matrix].any? { |c| other.is_a? c } other else ulysses20200426121236_add other end end end module Enumerable def sum init = 0, &block (block ? map(&block) : self).reduce init, :+ end end class Array def inner other zip(other).sum { |a, b| a b } end end (Again, note that this is Ruby 1.9.2.) Finally, just combine them up, and we can solve a Hamiltonian numerically. def solve_hamiltonian n, qp0, max_t, dt, hamiltonian runge_kutta qp0, max_t, dt, CLASSIC_4TH, ->t, qp do dqpdt = div qp, 1e-6, ->x { hamiltonian.(t, x) } Vector[dqpdt[n...n2], dqpdt[0...n].map(&:-@)] end end For example, let’s simulate a double pendulum. solve_hamiltonian 2,Vector[PI/2,0.0,0.0,0.0],Float::INFINITY,1e-3, ->t,(q1,q2,p1,p2){p12+p22/2+cos(q1-q2)p1p2-cos(q1)-cos(q2)} The codes are not complete in this post. See the attached file for details. You can open the project using RPG Maker VX Ace. The Game.exe file is not the official Game.exe executable but the third-party improved version of it called RGD (of version 1.3.2, while the latest till now is 1.5.1).
https://ulysseszh.github.io/physics/2020/04/27/simulation-rgss.html
Our goal is to simulate a mechanical system according to its Hamiltonian $\mathcal H\left(\mathbf q,\mathbf p,t\right)$. To utilize the canonical equations \begin{equation} \frac{\mathrm d\mathbf q}{\mathrm dt}= \frac{\partial\mathcal H}{\partial\mathbf p},\quad \frac{\mathrm d\mathbf p}{\mathrm dt}= -\frac{\partial\mathcal H}{\partial\mathbf p}, \label{canonical} \end{equation} we need to calculate the partial derivatives of $\mathcal H$. Here is a simple code to calculate partial derivatives. def div x0, dx, f f0 = f.(x0) n = x0.size Array.new n do |i| (f.(x0 + Vector.basis(n, i) dx) - f0) / dx end end (RGSS do not have matrix.rb, you can copy one from the attached file below.) Here x0 is a Vector, f is a call-able object as a function of vectors, dx is a small scalar which we are going to take 1e-6. Let x = Vector[q, p], and then Formula \ref{canonical} has the form \begin{equation} \frac{\mathrm d\mathbf x}{\mathrm dt}=f\left(\mathbf x\right). \end{equation} To solve this equation numerically, we need to use a famous method called the (explicit) Runge–Kutta method. def runge_kutta initial, max_t, dt, (pyramid, coefs), func (0..max_t).step(dt).reduce initial do |ret, t| $canvas.trace t, ret if $canvas coefs.zip(pyramid).each_with_object([]).sum do |(coef, row), ary| coef ary.push(func.(t, row.inner(ary) dt + ret)).last end dt + ret#p(ret) end end Note that Runge–Kutta is a family of methods. The argument (pyramid, coefs) takes one of the following, each of which is a single Runge–Kutta method. FORWARD_EULER = [[],[1]] EXPLICIT_MIDPOINT = [[],[1/2.0],[0,1]] HEUN = [[],[1],[1/2.0,1/2.0]] RALSTON = [[],[2/3.0],[1/4.0,3/4.0]] KUTTA_3RD = [[],[1/2.0],[-1,2],[1/6.0,2/3.0,1/6.0]] HEUN_3RD = [[],[1/3.0],[0,2/3.0],[1/4.0,0,3/4.0]] RALSTON_3RD = [[],[1/2.0],[0,3/4.0],[2/9.0,1/3.0,4/9.0]] SSPRK3 = [[],[1],[1/4.0,1/4.0],[1/6.0,1/6.0,2/3.0]] CLASSIC_4TH = [[],[1/2.0],[0,1/2.0],[0,0,1],[1/6.0,1/3.0,1/3.0,1/6.0]] RALSTON_4TH = [[],[0.4],[0.29697761,0.15875964],[0.21810040,-3.05096516, 3.83286476],[0.17476028, -0.55148066, 1.20553560, 0.17118478]] THREE_EIGHTH_4TH = [[],[1/3.0],[-1/3.0,1],[1,-1,1],[1/8.0,3/8.0,3/8.0,1/8.0]] Here we are going to take CLASSIC_4TH. The $canvas appearing here is an object that is going to draw the result onto the screen. Here we also need to have some patches to get it work. class Float alias ulysses20200426121236_add + def + other if zero? && [Vector, Matrix].any? { |c| other.is_a? c } other else ulysses20200426121236_add other end end end module Enumerable def sum init = 0, &block (block ? map(&block) : self).reduce init, :+ end end class Array def inner other zip(other).sum { |a, b| a b } end end (Again, note that this is Ruby 1.9.2.) Finally, just combine them up, and we can solve a Hamiltonian numerically. def solve_hamiltonian n, qp0, max_t, dt, hamiltonian runge_kutta qp0, max_t, dt, CLASSIC_4TH, ->t, qp do dqpdt = div qp, 1e-6, ->x { hamiltonian.(t, x) } Vector[dqpdt[n...n2], dqpdt[0...n].map(&:-@)] end end For example, let’s simulate a double pendulum. solve_hamiltonian 2,Vector[PI/2,0.0,0.0,0.0],Float::INFINITY,1e-3, ->t,(q1,q2,p1,p2){p12+p22/2+cos(q1-q2)p1p2-cos(q1)-cos(q2)} The codes are not complete in this post. See the attached file for details. You can open the project using RPG Maker VX Ace. The Game.exe file is not the official Game.exe executable but the third-party improved version of it called RGD (of version 1.3.2, while the latest till now is 1.5.1).