用 Ruby 做数值模拟和可视化

有丘直方

昨天下午的产物：
[pre lang="ruby"]require 'matrix'
require 'bigdecimal'
require 'openrgss'

RGSS.title = 'Runge-Kutta'
RGSS.show_fps = true

module NumberRefinements
        Float.prepend self
        Integer.prepend self
        BigDecimal.prepend self
        def + other
                zero? && [Vector, Matrix].any? { other.is_a? _1 } ? other : super
        end
end

class Vector
        def self.build(...)
                new Array.new(...)
        end
        def to_ary
                to_a
        end
end

class Array
        def inner other
                zip(other).sum { _1 * _2 }
        end
end

METHODS = [
FORWARD_EULER = [[],[1]],
EXPLICIT_MIDPOINT = [[],[1/2r],[0,1]],
HEUN = [[],[1],[1/2r,1/2r]],
RALSTON = [[],[2/3r],[1/4r,3/4r]],
KUTTA_3RD = [[],[1/2r],[-1,2],[1/6r,2/3r,1/6r]],
HEUN_3RD = [[],[1/3r],[0,2/3r],[1/4r,0,3/4r]],
RALSTON_3RD = [[],[1/2r],[0,3/4r],[2/9r,1/3r,4/9r]],
SSPRK3 = [[],[1],[1/4r,1/4r],[1/6r,1/6r,2/3r]],
CLASSIC_4TH = [[],[1/2r],[0,1/2r],[0,0,1],[1/6r,1/3r,1/3r,1/6r]],
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/3r],[-1/3r,1],[1,-1,1],[1/8r,3/8r,3/8r,1/8r]]
]

class Canvas
        attr_accessor :bitmap, :color
        def initialize
                @bitmap = Bitmap.new 1024, 768
                @color = Color.new 255, 255, 255
        end
        def trace t, ret
                block = ->y { @bitmap.set_pixel t*128, 200*(y+2), @color }
                ret.is_a?(Enumerable) ? ret.each(&block) : block.(ret)
        end
end

$canvas = Canvas.new

def runge_kutta initial, max_t, dt, (*pyramid, coefs), func
        (0..max_t).step(dt).reduce initial do |ret, t|
                $canvas&.trace t, ret
                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

def div x0, dx, f
        f0 = f.(x0)
        n = x0.size
        Vector.build n do |i|
                (f.(x0 + Vector.basis(size: n, index: i) * dx) - f0) / dx
        end
end

def solve_hamiltonian n, x0, t, dt, hamiltonian
        runge_kutta Vector[*x0], t, dt, CLASSIC_4TH, ->t, x do
                dxdt = div x, 1e-5, ->qp { hamiltonian.(t, qp) }
                Vector[*dxdt[n...n*2], *dxdt[0...n].map(&:-@)]
        end
end

BigDecimal.limit 25

include RGSS

Graphics.resize_screen 1024, 768
rgss_main do
        Thread.start do
                Sprite.new.bitmap = $canvas.bitmap
                rgss_stop
        end
end

alias ` BigDecimal

#METHODS.each do |method|
#        p runge_kutta `1`, `1`, `1e-4`, method, ->t,x { x*t }
#end

#p runge_kutta Vector[`1`,`0`], `4`, `1e-4`, CLASSIC_4TH, ->t,v { v.cross }
p solve_hamiltonian 1, Vector[`0`,`1`], `8`, `1e-4`, ->t,(x,p) { x**2+0.3*x**3 + p**2 }
gets
[/pre]
这个 openrgss 是我魔改过的：https://github.com/UlyssesZh/openrgss
这里模拟了个非简谐振动，Hamiltonian 是 H=x^2+0.3x^3+p^2
看下效果：

这里有个 RGSS 版的：
RungeKutta.rar (1.25 MB, 下载次数: 14)

2020-4-26 13:07 上传

点击文件名下载附件

---

把 RGSS 版的稍微修改了一下。这里做了个双摆的模拟。

RungeKutta.rar (1.25 MB, 下载次数: 27)

2020-4-27 13:34 上传

点击文件名下载附件