Project1

标题: 用 Ruby 做数值模拟和可视化 [打印本页]
作者: 有丘直方    时间: 2020-4-26 13:11
标题: 用 Ruby 做数值模拟和可视化
本帖最后由 有丘直方 于 2020-4-27 13:34 编辑

昨天下午的产物：
RUBY 代码复制
  1. require 'matrix'
  2. require 'bigdecimal'
  3. require 'openrgss'
  4.  
  5. RGSS.title = 'Runge-Kutta'
  6. RGSS.show_fps = true
  7.  
  8. module NumberRefinements
  9.         Float.prepend self
  10.         Integer.prepend self
  11.         BigDecimal.prepend self
  12.         def + other
  13.                 zero? && [Vector, Matrix].any? { other.is_a? _1 } ? other : super
  14.         end
  15. end
  16.  
  17. class Vector
  18.         def self.build(...)
  19.                 new Array.new(...)
  20.         end
  21.         def to_ary
  22.                 to_a
  23.         end
  24. end
  25.  
  26. class Array
  27.         def inner other
  28.                 zip(other).sum { _1 * _2 }
  29.         end
  30. end
  31.  
  32. METHODS = [
  33. FORWARD_EULER = [[],[1]],
  34. EXPLICIT_MIDPOINT = [[],[1/2r],[0,1]],
  35. HEUN = [[],[1],[1/2r,1/2r]],
  36. RALSTON = [[],[2/3r],[1/4r,3/4r]],
  37. KUTTA_3RD = [[],[1/2r],[-1,2],[1/6r,2/3r,1/6r]],
  38. HEUN_3RD = [[],[1/3r],[0,2/3r],[1/4r,0,3/4r]],
  39. RALSTON_3RD = [[],[1/2r],[0,3/4r],[2/9r,1/3r,4/9r]],
  40. SSPRK3 = [[],[1],[1/4r,1/4r],[1/6r,1/6r,2/3r]],
  41. CLASSIC_4TH = [[],[1/2r],[0,1/2r],[0,0,1],[1/6r,1/3r,1/3r,1/6r]],
  42. RALSTON_4TH = [[],[0.4],[0.29697761,0.15875964],[0.21810040,-3.05096516,3.83286476],[0.17476028, -0.55148066, 1.20553560, 0.17118478]],
  43. THREE_EIGHTH_4TH = [[],[1/3r],[-1/3r,1],[1,-1,1],[1/8r,3/8r,3/8r,1/8r]]
  44. ]
  45.  
  46. class Canvas
  47.         attr_accessor :bitmap, :color
  48.         def initialize
  49.                 @bitmap = Bitmap.new 1024, 768
  50.                 @color = Color.new 255, 255, 255
  51.         end
  52.         def trace t, ret
  53.                 block = ->y { @bitmap.set_pixel t*128, 200*(y+2), @color }
  54.                 ret.is_a?(Enumerable) ? ret.each(&block) : block.(ret)
  55.         end
  56. end
  57.  
  58. $canvas = Canvas.new
  59.  
  60. def runge_kutta initial, max_t, dt, (*pyramid, coefs), func
  61.         (0..max_t).step(dt).reduce initial do |ret, t|
  62.                 $canvas&.trace t, ret
  63.                 coefs.zip(pyramid).each_with_object([]).sum do |(coef, row), ary|
  64.                         coef * ary.push(func.(t, row.inner(ary) * dt + ret)).last
  65.                 end * dt + ret#p(ret)
  66.         end
  67. end
  68.  
  69. def div x0, dx, f
  70.         f0 = f.(x0)
  71.         n = x0.size
  72.         Vector.build n do |i|
  73.                 (f.(x0 + Vector.basis(size: n, index: i) * dx) - f0) / dx
  74.         end
  75. end
  76.  
  77. def solve_hamiltonian n, x0, t, dt, hamiltonian
  78.         runge_kutta Vector[*x0], t, dt, CLASSIC_4TH, ->t, x do
  79.                 dxdt = div x, 1e-5, ->qp { hamiltonian.(t, qp) }
  80.                 Vector[*dxdt[n...n*2], *dxdt[0...n].map(&:-@)]
  81.         end
  82. end
  83.  
  84. BigDecimal.limit 25
  85.  
  86. include RGSS
  87.  
  88. Graphics.resize_screen 1024, 768
  89. rgss_main do
  90.         Thread.start do
  91.                 Sprite.new.bitmap = $canvas.bitmap
  92.                 rgss_stop
  93.         end
  94. end
  95.  
  96. alias ` BigDecimal
  97.  
  98. #METHODS.each do |method|
  99. #        p runge_kutta `1`, `1`, `1e-4`, method, ->t,x { x*t }
  100. #end
  101.  
  102. #p runge_kutta Vector[`1`,`0`], `4`, `1e-4`, CLASSIC_4TH, ->t,v { v.cross }
  103. p solve_hamiltonian 1, Vector[`0`,`1`], `8`, `1e-4`, ->t,(x,p) { x**2+0.3*x**3 + p**2 }
  104. gets

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

这里有个 RGSS 版的：
RungeKutta.rar (1.25 MB, 下载次数: 14)
[line]10[/line]
把 RGSS 版的稍微修改了一下。这里做了个双摆的模拟。

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



欢迎光临 Project1 (https://rpg.blue/) Powered by Discuz! X3.1