Ruby script to draw a lorenz attractor with Runge-Kutta's method.

lorenz_attractor_runge_kutta.rb

#! /usr/local/bin/ruby
# *************************************
# Lorenz attractor (Runge-Kutta method)
# *************************************
#
require 'numo/gnuplot'

class LorenzAttractorRungeKutta
  DT            = 1e-3     # Differential interval
  STEP          = 100000   # Time step count
  X_0, Y_0, Z_0 = 1, 1, 1  # Initial values of x, y, z

  def initialize
    @res = [[], [], []]
  end

  def exec
    xyz = [X_0, Y_0, Z_0]
    STEP.times do
      k_0 = lorenz(xyz)
      k_1 = lorenz(xyz.zip(k_0).map { |x, k| x + k * DT / 2.0 })
      k_2 = lorenz(xyz.zip(k_1).map { |x, k| x + k * DT / 2.0 })
      k_3 = lorenz(xyz.zip(k_2).map { |x, k| x + k * DT })
      3.times do |i|
        xyz[i] += (k_0[i] + 2 * k_1[i] + 2 * k_2[i] + k_3[i]) * DT / 6.0
        @res[i] << xyz[i]
      end
    end
    plot
  rescue => e
    $stderr.puts "[#{e.class}] #{e.message}"
    e.backtrace.each { |tr| $stderr.puts "\t#{tr}" }
    exit 1
  end

  private

  def lorenz(xyz, p=10, r=28, b=8/3.0)
    return [
      -p * xyz[0] + p * xyz[1],
      -xyz[0] * xyz[2] + r * xyz[0] - xyz[1],
      xyz[0] * xyz[1] - b * xyz[2]
    ]
  rescue => e
    raise
  end

  def plot
    x, y, z = @res

    begin
      Numo.gnuplot do
        set   terminal: "png"
        set   output:   "lorenz_attractor_runge_kutta.png"
        set   title:    "Lorenz attractor (Runge-Kutta method)"
        set   xlabel:   "x"
        set   ylabel:   "y"
        set   zlabel:   "z"
        unset :key
        splot x, y, z, with: :lines, linecolor: {rgb: "blue"}
      end
    rescue => e
      raise
    end
  end
end

LorenzAttractorRungeKutta.new.exec if __FILE__ == $0