Skip to content

Instantly share code, notes, and snippets.

@komasaru

komasaru/spline_interpolation.rb

Last active April 19, 2018 06:01
Download ZIP

Revisions

  1. komasaru revised this gist Apr 19, 2018. 1 changed file with 13 additions and 13 deletions.
    26 changes: 13 additions & 13 deletions spline_interpolation.rb
    View file
    Original file line number Diff line number Diff line change
    @@ -145,17 +145,17 @@ def plot
    end
    end

    exit 0 unless __FILE__ = $0

    # グラフ用配列
    x_g, y_g = Array.new, Array.new
    # ３次スプライン補間
    obj = SplineInterpolation.new
    X[0].step(X[-1], 0.1) do |x|
    y = obj.interpolate(x)
    puts "%8.4f, %8.4f" % [x, y]
    x_g << x
    y_g << y
    if __FILE__ == $0
    # グラフ用配列
    x_g, y_g = Array.new, Array.new
    # ３次スプライン補間
    obj = SplineInterpolation.new
    X[0].step(X[-1], 0.1) do |x|
    y = obj.interpolate(x)
    puts "%8.4f, %8.4f" % [x, y]
    x_g << x
    y_g << y
    end
    # グラフ描画
    Graph.new(x_g, y_g).plot
    end
    # グラフ描画
    Graph.new(x_g, y_g).plot
  2. komasaru revised this gist Jan 9, 2018. 1 changed file with 23 additions and 39 deletions.
    62 changes: 23 additions & 39 deletions spline_interpolation.rb
    View file
    Original file line number Diff line number Diff line change
    @@ -1,5 +1,4 @@
    #!/usr/local/bin/ruby
    # coding: utf-8
    #*********************************************
    # ３次スプライン補間
    # （gnuplot によるグラフ描画付き）
    @@ -27,11 +26,10 @@ def initialize

    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
    return @a[i] * (t - @x[i]) ** 3 \
    + @b[i] * (t - @x[i]) ** 2 \
    + @c[i] * (t - @x[i]) \
    + @d[i]
    rescue => e
    $stderr.puts "[#{e.class}] #{e.message}"
    exit 1
    @@ -40,28 +38,23 @@ def interpolate(t)
    private

    def calc_h
    h = Array.new
    0.upto(@n - 1) { |i| h << @x[i + 1] - @x[i] }
    return h
    return (0..@n - 1).map { |i| @x[i + 1] - @x[i] }
    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])
    return (1..@n -1).map do |i|
    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]
    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
    @@ -71,45 +64,38 @@ def gauss_jordan(matrix)
    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 }
    k.upto(n) { |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
    next if i == k
    d = matrix[i][k]
    k.upto(n) { |j| matrix[i][j] -= d * matrix[k][j] }
    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
    return (0..@n - 1).map { |i| (v[i + 1] - v[i]) / (6 * (@x[i + 1] - @x[i])) }
    end

    def calc_b(v)
    b = Array.new
    0.upto(@n - 1) { |i| b << v[i] / 2.0 }
    return b
    return (0..@n - 1).map { |i| v[i] / 2.0 }
    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
    return (0..@n - 1).map do |i|
    (@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]
    return @y
    end

    def search_i(t)
    i, j = 0, j = @x.size - 1
    i, j = 0, @x.size - 1
    while i < j
    k = (i + j) / 2
    if @x[k] < t
    @@ -137,15 +123,13 @@ def plot
    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"
    @@ -161,17 +145,17 @@ def plot
    end
    end

    exit 0 unless __FILE__ = $0

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

    # ３次スプライン補間
    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)
    puts "%8.4f, %8.4f" % [x, y]
    x_g << x
    y_g << y
    end

    # グラフ描画
    Graph.new(x_g, y_g).plot
  3. komasaru created this gist Dec 7, 2015.
    177 changes: 177 additions & 0 deletions spline_interpolation.rb
    View file
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,177 @@
    #!/usr/local/bin/ruby
    # coding: utf-8
    #*********************************************
    # ３次スプライン補間
    # （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

    # ３次スプライン補間
    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