mk-mode BLOG

このブログは自作の自宅サーバに構築した Debian GNU/Linux で運用しています。
PC・サーバ構築等の話題を中心に公開しております。(クローンサイト: GitHub Pages

ブログ開設日2009-01-05
サーバ連続稼働時間
Reading...
Page View 合計
Reading...
今日
Reading...
昨日
Reading...

Ruby - 3次スプライン補間!

[ プログラミング, 数学 ] [ Ruby ]

こんばんは。

過去に「ラグランジュ補間」や「ニュートン補間」を C++ や Ruby で実装したことがありました。

今回は「3次スプライン補間」を Ruby で実装してみました。

0. 前提条件

  • Ruby 2.2.3-p173 での作業を想定。
  • グラフも描画するので、RubyGems ライブラリの gnuplot をインストールしておく。

1. 3次スプライン補間について

こちらの内容を自分なりに理解して以下のようにまとめた。

(数式が多いので別途 で記載した文書を貼り付け)

SPLINE_INTERPOLATION_1 SPLINE_INTERPOLATION_2 SPLINE_INTERPOLATION_3 SPLINE_INTERPOLATION_4 SPLINE_INTERPOLATION_5 SPLINE_INTERPOLATION_6

(上記文書の PDF はこちら

2. Ruby スクリプトの作成

以下のような Ruby スクリプトを作成する。
ロジックは前項の説明(アルゴリズム)どおりなので説明しない。(ただ、連立一次方程式の解法には「ガウス・ジョルダン法」を使用した)

spline_interpolation.rb
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
#!/usr/local/bin/ruby
# coding: utf-8
#*********************************************
# 3次スプライン補間
# (gnuplot によるグラフ描画付き)
#*********************************************
#
require 'gnuplot'

# N + 1 個の点
X = [0.0, 2.0, 3.0, 5.0, 7.0, 8.0]
Y = [0.8, 3.2, 2.8, 4.5, 2.5, 1.9]

class SplineInterpolation
  def initialize
    @x, @y = X, Y
    @n = @x.size - 1
    h = calc_h
    w = calc_w(h)
    matrix = gen_matrix(h, w)
    v = [0, gauss_jordan(matrix), 0].flatten
    @b = calc_b(v)
    @a = calc_a(v)
    @d = calc_d(v)
    @c = calc_c(v)
  end

  def interpolate(t)
    i = search_i(t)
    s  = @a[i] * (t - @x[i]) ** 3
    s += @b[i] * (t - @x[i]) ** 2
    s += @c[i] * (t - @x[i])
    s += @d[i]
    return s
  rescue => e
    $stderr.puts "[#{e.class}] #{e.message}"
    exit 1
  end

  private

  def calc_h
    h = Array.new
    0.upto(@n - 1) { |i| h << @x[i + 1] - @x[i] }
    return h
  end

  def calc_w(h)
    w = Array.new
    1.upto(@n - 1) do |i|
      w << 6 * ((@y[i + 1] - @y[i]) / h[i] - (@y[i] - @y[i - 1]) / h[i - 1])
    end
    return w
  end

  def gen_matrix(h, w)
    matrix = Array.new(@n - 1).map { Array.new(@n, 0) }
    0.upto(@n - 2) do |i|
      matrix[i][i]     = 2 * (h[i] + h[i + 1])
      matrix[i][-1]    = w[i]
      next if i == 0
      matrix[i - 1][i] = h[i]
      matrix[i][i - 1] = h[i]
      matrix[i][i - 1] = h[i]
    end
    return matrix
  end

  def gauss_jordan(matrix)
    v = Array.new
    n = @n - 1
    0.upto(n - 1) do |k|
      p = matrix[k][k]
      k.upto(n) { |j| matrix[k][j] = matrix[k][j] / p.to_f }
      0.upto(n - 1) do |i|
        unless i == k
          d = matrix[i][k]
          k.upto(n) { |j| matrix[i][j] = matrix[i][j] - d * matrix[k][j] }
        end
      end
    end
    matrix.each { |row| v << row[-1] }
    return v
  end

  def calc_a(v)
    a = Array.new
    0.upto(@n - 1) { |i| a << (v[i + 1] - v[i]) / (6 * (@x[i + 1] - @x[i])) }
    return a
  end

  def calc_b(v)
    b = Array.new
    0.upto(@n - 1) { |i| b << v[i] / 2.0 }
    return b
  end

  def calc_c(v)
    c = Array.new
    0.upto(@n - 1) do |i|
      c << (@y[i + 1] - @y[i]) / (@x[i + 1] - @x[i]) \
         - (@x[i + 1] - @x[i]) * (2 * v[i] + v[i + 1]) / 6
    end
    return c
  end

  def calc_d(v)
    return @y[0..-1]
  end

  def search_i(t)
    i, j = 0, j = @x.size - 1
    while i < j
      k = (i + j) / 2
      if @x[k] < t
        i = k + 1
      else
        j = k
      end
    end
    i -= 1 if i > 0
    return i
  end
end

class Graph
  def initialize(x_g, y_g)
    @x_g, @y_g, @x, @y = x_g, y_g, X, Y
  end

  def plot
    Gnuplot.open do |gp|
      Gnuplot::Plot.new(gp) do |plot|
        plot.terminal "png enhanced font 'IPA P ゴシック' fontscale 1.2"
        plot.output   "spline_interpolation.png"
        plot.title    "スプライン補間"
        plot.xlabel   "x"
        plot.ylabel   "y"
        plot.grid

        # 計算によって得られた点
        plot.data << Gnuplot::DataSet.new([@x_g, @y_g]) do |ds|
          ds.with      = "points"
          ds.linewidth = 2
          ds.linecolor = 3
          ds.notitle
        end

        # 予め与えらた点
        plot.data << Gnuplot::DataSet.new([@x, @y]) do |ds|
          ds.with      = "points"
          ds.linewidth = 2
          ds.linecolor = 1
          ds.notitle
        end
      end
    end
  rescue => e
    $stderr.puts "[#{e.class}] #{e.message}"
    exit 1
  end
end

# グラフ用配列
x_g, y_g = Array.new, Array.new

# 3次スプライン補間
obj = SplineInterpolation.new
X[0].step(X[-1], 0.1) do |x|
  y = obj.interpolate(x)
  printf("%8.4f, %8.4f\r\n", x, y)
  x_g << x
  y_g << y
end

# グラフ描画
Graph.new(x_g, y_g).plot

3. Ruby スクリプトの実行

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
$ ./spline_interpolation.rb
  0.0000,   0.8000
  0.1000,   0.9861
  0.2000,   1.1713
  0.3000,   1.3544
  0.4000,   1.5346
  0.5000,   1.7108
  0.6000,   1.8820
  0.7000,   2.0473
  0.8000,   2.2056
  0.9000,   2.3559
  1.0000,   2.4973
  1.1000,   2.6287
  1.2000,   2.7492
  1.3000,   2.8578
  1.4000,   2.9534
  1.5000,   3.0351
  1.6000,   3.1019
  1.7000,   3.1528
  1.8000,   3.1868
  1.9000,   3.2028
  2.0000,   3.2000
  2.1000,   3.1782
     :         :
====< 途中省略>====
     :         :
  7.5000,   2.1279
  7.6000,   2.0754
  7.7000,   2.0275
  7.8000,   1.9831
  7.9000,   1.9410
  8.0000,   1.9000

4. グラフ確認

Ruby スクリプトと同じディレクトリ内に “spline_interpolation.png” という画像ファイルが存在するはずなので、確認してみる。
(赤色の x が予め与えられた点、水色の + が補間された点)

SPLINE_INTERPOLATION

5. 参考サイト


スプライン補間ができないと先に進めない(私的な)事案があったために今回まとめてみた次第です。

以上。

Comments