Loading [MathJax]/jax/output/CommonHTML/jax.js

スプライン曲線の計算方法

2024年01月15日(月) 13時05分更新 icon 項目のみ表示/展開表示の切り替え

3次スプライン曲線の作成方法

既知点

下図のように6点の既知点からなる3次スプライン曲線を例に算出方法を説明する。
既知点をxn,ynとすると座標一覧は下表のとおりである。
XY
00 2041.68 1575.59
1 1 2830.84 2298.99
2 2 3685.07 1900.71
3 3 4191.15 1353.67
N-1 4 4889.14 1648.74
N 5 5403.39 2536.84

3次スプライン曲線による補間の基本

3次スプライン曲線による補間は、各区間(例えば(x 0,y 0)～ (x 1,y 1)区間など)事に3次方程式の係数 an , bn , cn , dnを算出して描画することである。
三次方程式は下式のとおりである。
(xn,yn)～(xn+1),(yn+1)区間 既知点が0からNまでの場合、式はN-1の補間式が作成できる。
今回の6点が既知点の場合
0 (x0,y0)～(x1),(y1)区間S0=a0(x−x0)3+b0(x−x0)2+c0(x−x0)+d0
1 (x1,y1)～(x2),(y2)区間S1=a1(x−x1)3+b1(x−x1)2+c1(x−x1)+d1
2 (x2,y2)～(x3),(y3)区間S2=a2(x−x2)3+b2(x−x2)2+c2(x−x2)+d2
3 (x3,y3)～(x4),(y4)区間S3=a3(x−x3)3+b3(x−x3)2+c3(x−x3)+d3
4 (x4,y4)～(x5),(y5)区間S4=a4(x−x4)3+b4(x−x4)2+c4(x−x4)+d4

各係数の算定

dnの算定
各区間のスプライン曲線の始終点はすべて既知点を通る。したがって、区間の補間式の両端のxに対するyの値は既知である。各区間の始点側はx=xnなので(x-xn)=0となる。よって yn=Sn=an(x−xn)3+bn(x−xn)2+cn(x−xn)+dn
dn=yn 式-1
bnの算定
隣同士の区間の曲線がスムーズに接続されるには、隣同士の曲線傾きが同一である必要がある。傾きは微分値なので 隣同士の区間の一次導関数、二次導関数の値が同一である必要がある。 とりあえずSn(x)に対する一次・二次導関数を求める Sn=an(x−xn)3+bn(x−xn)2+cn(x−xn)+dn
べき乗img alt="f(x)=x^n" class="img2latex">に対する一次導関数はf′(x)=n⋅x(n−1)なので、
一次導関数S′n(x−xn)=3an(x−xn)2+2bn(x−xn)+cn
二次導関数は、一次導関数をさらに微分したものなので、
S″(x−xn)=6an(x−xn)+2bn
始点側はx=xnなので(x-xn)=0となる。よって
S″(x−xn)=6an(x−xn)+2bn
bn=S″(x−xn)2 式-2

anの算定

隣同士の二次導関数の値は等しいので、
S″(x−xn)=6an(x−xn)+2bn
an=(S″n+1(x)−2bn)6(xn+1−xn)
an=(S″n+1(x)−S″n(x))6(xn+1−xn)式-3
cnの算定
補間式は、与点を必ず通るので、
an(x−xn)3+bn(x−xn)2+cn(x−xn)+dn=yn+1
式-1,2,3を代入すると
S″n+1(x)−2bn6(xn+1−xn)(xn+1−xn)(xn+1−xn)3+S″n(x)2(xn+1−xn)2+Cn(xn+1−xn)=Yn=Yn+1
cnについて解くと
cn(xn+1−xn)=Yn+1−Yn−S″n+1(x)−S″n6(xn+1−xn)(xn+1−xn)(xn+1−xn)3+S″n(x)2(xn+1−xn)2
cn=Yn+1−Yn−S″n+1(x)−S″n6(xn+1−xn)(xn+1−xn)(xn+1−xn)3+S″n(x)2(xn+1−xn)2xn+1−xn
cn=Yn+1−Ynxn+1−xn−(S″n+1(x)−S″n(x))(xn+1−xn)26(xn+1−xn)−S″n(x)(xn+1−xn)2
上記式の一部を整理するために下式の簡略化をする
−(S1−S)(x1−x)6−S(x1−x)6
展開すると
−S1x16−S˙x13+S1x16+S⋅x3
整理すると
−S1x1−2S⋅x1+2S1⋅x6
(a+b)(c−d)=−bd−ad+bc+acより
−16(x1−x)(2S+S1)
整理結果をもとの式に反映させる
cn=Yn+1−YnXn+1−xn−16(xn+1−xn)(2S″n(x)+2S″n+1(x))

二次導関数の算出

隣同士の一次導関数の値は一致する必要があるので
S′n(xn+1)=S″n+1(xn+1)
一次導関数は
S′n(x−xn+1)=3an(x−xn)2+2bn(x−xn)+cn
xn+1の時 x-xnは0となるので、
3a′n(xn+1−xn)2+2bn(xn−xn)+cn
上記式に an,bn,cnの式を代入する。
an=S″n+1(xn)−S″n(xn)6(xn+1−xn)
bn=S″n(xn)2
cn=(yn+1−yn)(xn+1−xn)−16(xn+1−xn)(2S″n(xn)+S″n+1(xn))
cn+1=(yn+2−yn+1)(xn+2−xn+1)−16(xn+2−xn+1)(2S″n(xn+1)+S″n+2(xn))
3(S″n+1(xn)−S″n(xn)6(xn+1−xn))(xn+1−xn)2+2(S″n(xn)2)(xn+1−xn)
+{(yn+1−yn)(xn+1−xn)−16(xn+1−xn)(2S″n(xn)+S″n+1(xn))}=
(yn+2−yn+1)(xn+2−xn+1)−16(xn+2−xn+1)(2S″n+1(xn+1)+S″n+2(xn))
Sn" Sn+1" Sn+2"を式の左側に移行する。
S″n+1(xn)−S″n(xn)2(xn+1−xn)+S″n(xn)(xn+1−xn)−16(xn+1−xn)(2S″n(xn)+S″n+1(xn))
+16(xn+2−xn+1)(2S″n+1(xn+1)+S″n+2(Xn))=(yn+2−yn+1)(xn+2−xn+1)−(yn+1−yn)(xn+1−xn)
式を分解する
S″n+1(xn)2(xn+1−xn)−S″n(xn)(xn+1−xn)−13(xn+1−xn)S″n(xn)
−16(xn+1−xn)S″n+1(xn)+13(xn+2−xn+1)S″n+1(xn+1)
+16(xn+2−xn+1)S″n+2(xn)=(yn+2−yn+1(xn+2−xn+1)−(yn+1−yn)(xn+1−xn)
式を整理する
S″n(xn){−(xn+1−xn)2+(xn+1−xn)−13(xn+1−xn)}
S″n+1(xn){−12(xn+1−xn)−16(xn+1−xn)−13(xn+2−xn+1)}+16(xn+2−xn+1)S″n+2(xn)=
(yn+2−yn+1)(xn+2−xn+1)−(yn+1−yn)(xn+1−xn)
S″n(xn)(xn+1−xn)6+S″n+1(xn){13(xn+2−xn+1)}
+16(xn+2−xn+1)S″n+2(xn)=(yn+2−yn+1)xn+2−xn+1−(yn+1−yn)(xn+1−xn)
式を6倍する
(xn+1−xn)Sn″(xn)+2(xn+2−xn+1)S″n+1(xn)+(xn+2−xn+1)S″n+2(xn)=6{(yn+2−yn+1)(xn+2−xn+1)−(yn+1−yn)(xn+1−xn)}
各区間の二次導関数について連立方程式を作成する。
今回の6点が既知点の場. 0 (x0,y0)～(x1),(y1)区間
n=0
(xn+1−xn)Sn″(xn)+2(xn+2−xn+1)S″n+1(xn)+(xn+2−xn+1)S″n+2(xn)=6{(yn+2−yn+1)(xn+2−xn+1)−(yn+1−yn)(xn+1−xn)}
1 (x1,y1)～(x2),(y2)区間
n=1
(xn+1−xn)Sn″(xn)+2(xn+2−xn+1)S″n+1(xn)+(xn+2−xn+1)S″n+2(xn)=6{(yn+2−yn+1)(xn+2−xn+1)−(yn+1−yn)(xn+1−xn)}
2 (x2,y2)～(x3),(y3)区間
n=2
(xn+1−xn)Sn″(xn)+2(xn+2−xn+1)S″n+1(xn)+(xn+2−xn+1)S″n+2(xn)=6{(yn+2−yn+1)(xn+2−xn+1)−(yn+1−yn)(xn+1−xn)}
3 (x3,y3)～(x4),(y4)区間
n=3
(xn+1−xn)Sn″(xn)+2(xn+2−xn+1)S″n+1(xn)+(xn+2−xn+1)S″n+2(xn)=6{(yn+2−yn+1)(xn+2−xn+1)−(yn+1−yn)(xn+1−xn)}
4 (x4,y4)～(x5),(y5)区間
n=4
(xn+1−xn)Sn″(xn)+2(xn+2−xn+1)S″n+1(xn)+(xn+2−xn+1)S″n+2(xn)=6{(yn+2−yn+1)(xn+2−xn+1)−(yn+1−yn)(xn+1−xn)}
一般的に表すと0～Nまでの既知点がある場合
0～N-2について下記式の連立方程式が成り立つ。
(xn+1−xn)Sn″(xn)+2(xn+2−xn+1)S″n+1(xn)+(xn+2−xn+1)S″n+2(xn)=6{(yn+2−yn+1)(xn+2−xn+1)−(yn+1−yn)(xn+1−xn)}
なお、S0"(x)=SN-1"(x)=0とする。（自然スプライン）
式を簡素化するために下式とおく。
hn=xn+1−xn
hn+1−hn=xn+2−xn+1+xn+1−xn=xn+2−xn
kn=6{yn+1−ynhn−yn−yn−1hn−1}
2(h1+h0)s″1+h1S″2=k1
h1S″1+2(h2+h1)S″2+h2S″3=K2
h2S″2+2(h3+h2)S″3=K2
より一般的な表記をすると0～N-2について下記式の連立方程式が成り立つ。
n=0の時 2(h1+h0)s″1+h1S″2=k1
n=1～N-4 hnS″n+2(hn+1+hn)S″n+1+hn+1S″n+2=Kn+1
n=N-3hN−3S″N−3+2(hN−2+hN−3)S″N−2=KN−3
例えば2つ変数がある連立方程式を行列表記すると下式となる。
3x−4y=−5
2x+y=4
[3−421][xy]=[−54]
解き方はいろいろあるがここではクラメールの公式の例とガウスの消去法を示す。
連立方式を行列で示すと下式のとおりとなる。
[2(h1+h0)h10h12(h2+h1)h20h22(h3+h2)] [S″1S″2S″3]=[k1k2k3]
　クラメール
[a11a12a21a22][xy]=[b1b2]
x=|b1a12b2a22||a11a12a21a22| y=|a11b1a21b2||a12a12a22a22|
なお、||は絶対値ではなく、行列の値を求めるという意味である。例えば 2次の正方行列の場.、下式で解を求めることができる。
|a11b1a21b2|=a11a22−a12a21
3次の正方行列の場.サラスの方法により求めることができる。 4次以上については各種方法がある。また、行列の値は、Excelおよびmaximax等で求めることができる。 Excelの場合、MDETERM関数を用いる S″1=|k1h10k22(h2+h1)h2k3h22(h3+h2)||2(h1+h0)h10h12(h2+h1)h20h22(h3+h2)|
S″2=|2(h1+h0)k10h1k2h20k32(h3+h2)||2(h1+h0)h10h12(h2+h1)h20h22(h3+h2)|
S″3=|2(h1+h0)h1k1h12(h2+h1)k20h2k3||2(h1+h0)h10h12(h2+h1)h20h22(h3+h2)|
各hn、knを算定する。
hn=xn+1−xn
kn=6(yn+1−ynhn−yn+1−ynhn)
nx y hn kn
0 2041.68 1575.59 789.16
1 2830.84 2298.99 854.23 -8.297492068
2 3685.07 1900.71 506.08 -3.688148198
3 4191.15 1353.67 697.99 9.0220696
4 4889.14 1648.74 514.25 7.825431565
5 5403.39 2536.84
[3286.78854.230854.232720.62506.080506.082408.14][S″1S″2S″3]=[−8.297492−3.6881489.0220696]
" 計算されたhn,knを代入する
S″1=|−8.297492068854.230−3.6881481982720.62506.089.0220696506.082408.14||3286.78854.230854.232720.62506.080506.082408.14|=−407497721.893×1010=−0.002152
S″2=|3286.78−8.2974920680854.23−3.688148198506.0809.02206962408.14||3286.78854.230854.232720.62506.080506.082408.14|=−271300471.893×1010=−0.001433
S″3=|3286.78854.23−8.2974921854.232720.62−3.68814820506.089.0220696||3286.78854.230854.232720.62506.080506.082408.14|=766402691.893×1010=0.0040476
ガウスの消去法
　　ガウスの前進消去法を用いて上三角行列を作成する。
[a11a12a13a21a22a23a31a32a33]
上記の行列２列目以降に下式を適用する。
ui1=ai1/a11
aji = aji - ui1 * a1j
bi=bi-ui1*b1
[a11a12a13a21−u11∗a11a22−u11∗a12a23−u11∗a13a31−u31∗a11a32−u31∗a12a33−u31∗a13]
式が複雑なので適用結果をaiｊ'と置く
[a11′a12′a13′0a22′a23′0a32′a33′]
更に3行目2列目以降の行に下式を適用する
ui2'=ai2'/a22'
aij'=aji'-ui2*a2j' bi=bi-ui1*bi
[a11′a12′a13′0a22′a23′0a32′=a32′−ui2′∗a22′a33′=a33′−ui2′∗a23′]
適用結果を整理すると下式となる
[a11′a12′a13′0a22′a23′00a33′=a33′−ui2′∗a23′]
以下、行がなくなるまで、順番に行と列を移動しながら適用していく。
上三角行列ができるので、左上から右下方向の対角線上の値を掛けていくと行列の値が算出できる。
a11'*a22'*a33'=a33'-ui2'*a23'
前進消去で行列式を上三角行列に変換する。
         1        2        3       kn
 1 3286.78   854.23        0  -8.2974921
 2  854.23  2720.62   506.08  -3.6881482
 3       0   506.08  2408.14   9.0220696
ui1=ai1/a11
aji = aji - ui1 * a1j
bi=bi-ui1*b1 
        1           2        3              ui1
1 3286.78      854.23        0  -8.2974921
2       0 2498.606694   506.08  -1.5316404  0.259898746
3       0      506.08  2408.14   9.0220696            0
aij'=aji'-ui2*a2j' bi=bi-ui1*bi
        1           2           3       kn     ui2'=ai2'/a22' 
1 3286.78      854.23           0  -8.2974921
2       0 2498.606694      506.08  -1.5316404
3       0           0 2305.636086  9.33229553  0.202544883
後退代入 S″3=9.3322955/2305.636086=0.004047601
S″2=(−1.53164+−506.08∗0.004047601)2498.606694=−0.001433
S″1=(=−8.297492+0∗0.004047601+−854.23∗−0.001432819)3286.78=−0.002152
S1"～S3"の値より an,bn,cn,dnを算定する。
an=S″n+1(x)−S″n(x))6(xn+1−xn)
bn=S″n(x))2
cn=(yn+1−yn)(xn+1−xn)−16(xn+1−xn)(2S″n(xn)+S″n+1(xn))
dn=yn
      x        y            Sn"           a           b             c         d
0  2041.68  1575.59             0  -4.54516E-07           0   1.199731674  1575.59
1  2830.84  2298.99  -0.002152117    1.4034E-07  -0.0010761   0.350549323  2298.99
2  3685.07  1900.71  -0.001432819   1.80486E-06  -0.0007164  -1.180630528  1900.71
3  4191.15  1353.67   0.004047601   -9.6649E-07   0.0020238  -0.518985955  1353.67
4  4889.14  1648.74             0             0           0    1.72698104  1648.74
5  5403.39  2536.84
xを変化させながら、an,bn,cn,dnより補間値ｙを算出する。
yn=Sn=an(x−xn)3+bn(x−xn)2+cn(x−xn)+dn
各区間のxを4分割して補間すると下表のとおりとなる。
     x         y          Sn"          a               b          c        d
     0  2041.68  1575.59
                                0  -4.54516E-07           0   1.199731674  1575.59
                                   -4.54516E-07           0   1.199731674  1575.59
                                   -4.54516E-07           0   1.199731674  1575.59
1  2830.84  2298.99  
                     -0.002152117    1.4034E-07  -0.0010761   0.350549323  2298.99
                                     1.4034E-07  -0.0010761   0.350549323  2298.99
                                     1.4034E-07  -0.0010761   0.350549323  2298.99
2  3685.07  1900.71
                     -0.001432819   1.80486E-06  -0.0007164  -1.180630528  1900.71
                                    1.80486E-06  -0.0007164  -1.180630528  1900.71
                                    1.80486E-06  -0.0007164  -1.180630528  1900.71
3  4191.15  1353.67
                      0.004047601   -9.6649E-07   0.0020238  -0.518985955  1353.67
                                    -9.6649E-07   0.0020238  -0.518985955  1353.67
                                    -9.6649E-07   0.0020238  -0.518985955  1353.67
4  4889.14  1648.74  
                                0             0           0    1.72698104  1648.74
                                              0           0    1.72698104  1648.74
                                              0           0    1.72698104  1648.74
5  5403.39  2536.84

      x-xn        x      x-xn        yn
0              2041.68                1575.59
               2238.97    197.29  1808.794746
               2436.26    394.58  2021.057593
               2633.55    591.87  2191.436644
1     789.16   2830.84                2298.99
             3044.3975  213.5575  2326.143715
              3257.955   427.115   2263.34747
             3471.5125  640.6725   2118.80249
2     854.23   3685.07                1900.71
               3811.59    126.52  1743.524121
               3938.11    253.04  1585.334364
               4064.63    379.56  1448.072425
3     506.08   4191.15                1353.67
             4365.6475  174.4975  1319.596429
              4540.145   348.995  1377.958061
             4714.6425  523.4925  1497.943163
4     697.99   4889.14                1648.74
             5017.7025  128.5625     1870.765
              5146.265   257.125      2092.79
             5274.8275  385.6875     2314.815
5     514.25   5403.39                2536.84
補間結果をグラフに表した結果が下図のとおりである。
PDF版