用 Ruby 做数值模拟和可视化

有丘直方 · 发表于 2020-4-26 13:11:09

加入我们，或者，欢迎回来。

您需要登录才可以下载或查看，没有帐号？注册会员

x

本帖最后由有丘直方于 2020-4-27 13:34 编辑

昨天下午的产物：

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

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

这个 openrgss 是我魔改过的：https://github.com/UlyssesZh/openrgss
这里模拟了个非简谐振动，Hamiltonian 是 H=x^2+0.3x^3+p^2
看下效果：

这里有个 RGSS 版的：

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

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

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

帐号		自动登录	找回密码
密码			注册会员

[原创发布] 用 Ruby 做数值模拟和可视化

加入我们，或者，欢迎回来。

评分