点をなめらかな曲線でつなぐ--3次式補間

一昨年のボーナスが80万円、今年のボーナスが75万円、去年のはいくらだったっけ？なんて時には大体その中間だろうと見当をつける。しかし、3年前に90万円だったりすると、去年はもう少し少なかったように思う。グラフを書いて見て、滑らかに線を引っ張ると半年前が底で去年も75万と言うことになる。

これをコンピュータでやろうとするとどうなるか？　一番簡単なのは2点の間を直線で結んで折れ線グラフにしてしまう線形補間だが、大体その中間だろうと見当をつけるのと同じで、これではあまりに芸がない。やはり点を滑らかにつなぎたいと言うことになる。

直線の次に簡単なのは放物線で、点が３つあれば1つの放物線が決められ、3点が綺麗に曲線で結ばれる。しかし、実際には3点でなく、もっとたくさん点がある場合が多い。この場合、区間に切って、区間ごとに放物線を決める3点を変えていくことになる。

やってみるとこれがうまく行かない。線がまがりくねっている場合、放物線を決める3点を変更したとたんに値にピョンと飛びが出てしまうのだ。点のところで切り替えをすれば飛びはなくなるが、カクっと角ができて、全然なめらかではない。放物線は2乗の係数の正負で曲がり方が逆になってしまうのだからいたし方ない。

世の中にはこれを解決する方法があり、ベジェ曲線とかスプライン曲線と言うものがそれだ。ベジェ曲線は2点間をあまりに自由につなぐので補間には向かない。スプライン曲線は数学的に完璧なものだが大きなマトリックスを解いたりせねばならずお手軽ではない。

実用的なのは放物線を一歩進めて3次曲線で結ぶという方法だ。補間する区間の端とその外側の点４つで3次式を決める。3次式の場合、4点をずらして行くと、区間の外側で傾きが反対になっても、区間内の曲がりは同じになる。微分連続というわけではないが、微分がそう急激に変わらずに連続性が保たれているわけだ。

(x(0),y(0))から(x(3),y(3))までの4点を通る3次曲線yy=ax^3+bx^2+cx+dの係数a,b,cを求めるのは以下のとおり。ちょっと複雑だが一度書いてしまえばいつでも滑らかに点をつなぐことができる。正式のスプラインは点列全部を読み込んで大きな行列を解くのだが、この3次式補間は周りの4点だけで値を決めるのではるかにお手軽になる。

強引に方程式を解いて ｙ=a+bx+cx^2+dx^3 の係数を決めてやると

d=(x(0)*(x(0)-x(2))*x(2)*(x(0)-x(3))*(x(2)-x(3))*x(3)*y(1)+x(1)^3*(x(0)*(x(0)-x(3))*x(3)*y(2)+x(2)^2*(x(3)*y(0)+x(0)*y(3))+x(2)*(x(3)^2*y(0)-x(0)^2*y(3)))+x(1)*(x(0)^2*(x(0)-x(3))*x(3)^2*y(2)+x(2)^3*(x(3)^2*y(0)+x(0)^2*y(3))+x(2)^2*(x(3)^3*y(0)-x(0)^3*y(3)))+x(1)^2*(x(0)*x(3)*(x(0)^2+x(3)^2)*y(2)+x(2)^3*(x(3)*y(0)-x(0)*y(3))+x(2)*(x(3)^3*y(0)+x(0)^3*y(3))))/((x(0)-x(1))*(x(0)-x(2))*(x(1)-x(2))*(x(0)-x(3))*(x(1)-x(3))*(x(2)-x(3)))
c=(x(0)^2*(x(0)-x(3))*x(3)^2*(y(1)-y(2))+x(2)^3*(x(3)^2*(y(0)-y(1))+x(0)^2*(y(1)-y(3)))+x(1)^2*(x(3)^3*(y(0)-y(2))+x(0)^3*(y(2)-y(3))+x(2)^3*(-y(0)+y(3)))+x(2)^2*(x(3)^3*(-y(0)+y(1))+x(0)^3*(-y(1)+y(3)))+x(1)^3*(x(3)^2*(-y(0)+y(2))+x(2)^2*(y(0)-y(3))+x(0)^2*(-y(2)+y(3))))/((x(0)-x(1))*(x(0)-x(2))*(x(1)-x(2))*(x(0)-x(3))*(x(1)-x(3))*(x(2)-x(3)))
b=(-x(0)*x(3)*(x(0)^2-x(3)^2)*(y(1)-y(2))+x(2)*(x(3)^3*(y(0)-y(1))+x(0)^3*(y(1)-y(3)))+x(1)^3*(x(3)*(y(0)-y(2))+x(0)*(y(2)-y(3))+x(2)*(-y(0)+y(3)))+x(2)^3*(x(3)*(-y(0)+y(1))+x(0)*(-y(1)+y(3)))+x(1)*(x(3)^3*(-y(0)+y(2))+x(2)^3*(y(0)-y(3))+x(0)^3*(-y(2)+y(3))))/((x(0)-x(1))*(x(0)-x(2))*(x(1)-x(2))*(x(0)-x(3))*(x(1)-x(3))*(x(2)-x(3)))
a=(x(0)*(x(0)-x(3))*x(3)*(y(1)-y(2))+x(2)^2*(x(3)*(y(0)-y(1))+x(0)*(y(1)-y(3)))+x(1)*(x(3)^2*(y(0)-y(2))+x(0)^2*(y(2)-y(3))+x(2)^2*(y(0)+y(3)))+x(2)*(x(3)^2*(y(0)+y(1))+x(0)^2*(y(1)+y(3)))+x(1)^2*(x(3)*(y(0)+y(2))+x(2)*(y(0)-y(3))+x(0)*(-y(2)+y(3))))/((x(0)-x(1))*(x(0)-x(2))*(x(1)-x(2))*(x(0)-x(3))*(x(1)-x(3))*(x(2)-x(3)))

なのだが、それは凡人がやることであって、賢く考えればもっと簡単に求めることができる。3次式による補間とはラグランジェ補間の次数を３次に止めただけのことだ。

y=y(4)*(x-x(1)/(x(4)-x(1)) * (x-x(2)/(x(4)-x(2)) * (x-x(3)/(x(4)-x(3))
+y(1)*(x-x(2)/(x(1)-x(2)) * (x-x(3)/(x(1)-x(3)) * (x-x(4)/(x(1)-x(4))
+y(2)*(x-x(2)/(x(2)-x(3)) * (x-x(4)/(x(2)-x(4)) * (x-x(1)/(x(2)-x(1))
+y(3)*(x-x(2)/(x(3)-x(4)) * (x-x(1)/(x(3)-x(1)) * (x-x(2)/(x(3)-x(2))

というのがその解なのだが、さすがラグランジェ先生は頭が良い。yの4つの値のミキシングだと考え、x=x(1)の時には他の項が消えるように(x-x(i))をかけて、あとで規格化するために割ってやる。これで何次の式であろうと作れる。ここでは最低必要な3次式にしている。これで4点ごとに3次式を決めて行けば補間ができて、緩やかな曲線なら本格的なスプラインに比べても遜色ないように思う。

「スプラインに比べて遜色ないように思う」と書いてはみたが、データの上がり下がりが大きい場合はやはり角張る。大きく外れるわけではないのだが厳密には微分連続ではないからだ。遜色ないようにするにはもう少し改良がいる。







[追記]

データの上がり下がりが大きい場合に角張るのはどうも気に入らない。大きく外れるわけではないのだが厳密には微分連続ではないからだ。遜色ないようにするにはもう少し改良がいる。必ずしも外側の2点を通らなくてもいいのだから、微分連続を確保したほうが良い。

「データ点での傾きはその前後の直線近時の傾きの平均とする」ということに決めてしまえば微分連続は区間が変わっても自動的に確保できる。この場合、

g(j) = y(j + 1) - y(j)
h(j) = x(j + 1) - x(j) として
a = y(i)
b = g(i - 1) / h(i - 1) / 2 + g(i) / h(i) / 2
c = -(g(1 + i) * h(i) - 5 * g(i) * h(1 + i) + 4 * b * h(i) * h(1 + i)) / (2 * h(i) ^ 2 * h(1 + i))
d = -(-g(1 + i) * h(i) + 3 * g(i) * h(1 + i) - 2 * b * h(i) * h(1 + i)) / (2 * h(i) ^ 3 * h(1 + i))

が係数となり、曲線の値は
yy = (((d * (xx - x(i)) + c) * (xx - x(i)) + b) * (xx - x(i)) + a)
で与えられる。この場合も4点だけで3次式が決められるので、エクセルの関数に組み込んで使うとお手軽だ。中の2点を通り、外の点は、通るのではなく、傾きを決めるのに使われていることになる。スプラインの場合二次微分まで連続とするのだが、そこまでせずとも微分連続だけで十分に滑らかといえる。