帰ってきたネコブログ

3次スプライン曲線の描き方について調べていたんだけど、
アルゴリズムと入力と出力の例みたいのを見つけたので書き直してみた。

参考にしたサイトはこちら。

３次スプライン補間法
http://next1.msi.sk.shibaura-it.ac.jp/MULTIMEDIA/numeanal/node16.html

use v5.14;
use strict;
use warnings;

my @points = (
    [ 0.9, 1.3  ],
    [ 1.3, 1.5  ],
    [ 1.9, 1.85 ],
    [ 2.1, 2.1  ],
    [ 2.6, 2.6  ],
    [ 3.0, 2.7  ],
    [ 3.9, 2.4  ],
    [ 4.4, 2.15 ],
    [ 4.7, 2.05 ],
    [ 5.0, 2.10 ],
    [ 6.0, 2.25 ],
    [ 7.0, 2.3  ],
    [ 8.0, 2.25 ],
    [ 9.2, 1.95 ]
);

my $n = scalar( @points ) - 1;
my @a_j = map { $_->[1]; } @points;

my @h_j = ();
for (my $i=0; $i<$n; $i++) {
    push @h_j, $points[$i+1]->[0] - $points[$i]->[0];
}

my @alpha_j = ( 0 );
for (my $i=1; $i<$n; $i++) {
    my $tmp  = ((3 / $h_j[$i]  ) * ($a_j[$i+1] - $a_j[$i  ]));
       $tmp -= ((3 / $h_j[$i-1]) * ($a_j[$i  ] - $a_j[$i-1]));
    push @alpha_j, $tmp;
}

my @l_j = ( 1 );
my @u_j = ( 0 );
my @z_j = ( 0 );

for (my $i=1; $i<$n; $i++) {
    push @l_j, ( (2 * ($points[$i+1]->[0] - $points[$i-1]->[0])) - ($h_j[$i-1] * $u_j[$i-1]) );
    push @u_j, $h_j[$i] / $l_j[$i];
    push @z_j, ($alpha_j[$i] - ($h_j[$i-1] * $z_j[$i-1])) / $l_j[$i];
}

push @l_j, ( 1 );
push @z_j, ( 0 );

my @b_j = map { 0; } 0..$n;
my @c_j = map { 0; } 0..$n;
my @d_j = map { 0; } 0..$n;
for (my $i=($n - 1); 0<=$i; $i--) {
    $c_j[$i] = $z_j[$i] - ($u_j[$i] * $c_j[$i+1]);
    $b_j[$i] = (($a_j[$i+1] - $a_j[$i]) / $h_j[$i]) - (($h_j[$i] / 3) * ($c_j[$i+1] + (2 * $c_j[$i])));
    $d_j[$i] = ($c_j[$i+1] - $c_j[$i]) / (3 * $h_j[$i]);
}

printf( "%2s, %6s, %6s, %6s, %6s, %6s\n",
    'i', 'x(i)', 'a(i)', 'b(i)', 'c(i)', 'd(i)' );
for (my $i=0; $i<=$n; $i++) {
    if ( $i == $n ) {
        printf( "%2d, %6.3f\n", $i, $points[$i]->[0] );
        next;
    }

    printf( "%2d, %6.3f, %6.3f, %6.3f, %6.3f, %6.3f\n",
        $i, $points[$i]->[0], $a_j[$i], $b_j[$i], $c_j[$i], $d_j[$i] );
}

$ perl aaa.pl 
 i,   x(i),   a(i),   b(i),   c(i),   d(i)
 0,  0.900,  1.300,  0.540,  0.000, -0.248
 1,  1.300,  1.500,  0.421, -0.297,  0.947
 2,  1.900,  1.850,  1.087,  1.407, -2.956
 3,  2.100,  2.100,  1.295, -0.367, -0.447
 4,  2.600,  2.600,  0.593, -1.037,  0.445
 5,  3.000,  2.700, -0.022, -0.502,  0.174
 6,  3.900,  2.400, -0.503, -0.032,  0.078
 7,  4.400,  2.150, -0.477,  0.085,  1.314
 8,  4.700,  2.050, -0.071,  1.268, -1.581
 9,  5.000,  2.100,  0.262, -0.155,  0.043
10,  6.000,  2.250,  0.080, -0.027, -0.003
11,  7.000,  2.300,  0.017, -0.036, -0.031
12,  8.000,  2.250, -0.147, -0.128,  0.036
13,  9.200

スプライン曲線のコードを読む時は、
データの数がN個なのか、それとも補間区間の数がN個なのか把握する必要があって、
検索して見つけたものの多くは、
補間区間をN個とし、データの数はN+1個だった。

今回、苦戦したのは、
変数名にa(A)とα(Alpha)が使われてて見分け付かなかった点で、
あと、後ろから係数を求めているのは「トーマス法」がキーワードっぽい。

とりあえず、おしまい。

カテゴリー: Perl, コンピューター by ネコ on 2015年03月27日 | タグ: スプライン補間