はじめに
実験結果をグラフに表示する時、測定点をプロットしてから、その間を自在定規などで結んだ経験があると思います。ここでは、これを自動的に行うプログラムを紹介します。このような作業は「補間」（Interpolation）と呼ばれますが、「補間」には色々な方式がありますので、これらを比較できるようにしました。
- 完成版アプレットを見る
対象読者
Lagrange、Newton、Splineなどの補間方式に興味を持ち、実験結果をまとめるのに、自分で作ったプログラムを応用したい人。また、CGの基礎としての補間の原理を学びたい人。
必要な環境
J2SE 5.0を使っていますが、それより古いバージョンでも大丈夫です。
補間とは
「狭義の補間」は、二点のデータを正確なものと仮定し、その中間点のデータを推測することで、「内挿」とも言います。これは、デジカメなどの画像処理で用いられます。一方、「広義の補間」は、与えられたすべてのデータ点を通るように曲線を描くこと、および、その曲線の方程式を求めることで、データ点に誤差があっても、そのまま補間されますので、補間は正しい値を推測するものではありません。補間は、実験結果のグラフ作成のほか、アウトラインフォント、CGなどに広く用いられています。
似たものに、「曲線の当てはめ（Curve Fitting）」があります。これは、あらかじめ与えられた次数の多項式で、測定点との誤差を最小にするもので、実験結果のグラフ化や実験式の作成に用いられます。この場合は、測定点を必ず通るのではなく、より正しい値を推測します。
補間のいろいろと、プログラムでの実現方法
0番目からn番目まで（n+1）個の点を必ず通る曲線を表現するには、n次の多項式が必要十分条件で、この多項式を求めることが補間です。
多項式を求める方法には、Lagrange多項式が良く知られていますが、Newtonの方法もあり、二つの方法の結果は同じになります。後者は、データ点の追加が容易な長所がある代わりに、再帰処理のために時間が掛かります。筆者の案によるLagrange補間の改良（すでに誰かがやっているかも知れませんが）も、併せて紹介します。
一方、多項式の次数を限定し、狭い区間ごとに異なる多項式を求める方法にSplineがあります。これは、曲線の円滑な連続性には欠けますが、異常な振動が起きません。Splineの中では、ここで紹介する3次Splineが代表的なものです。CGなどで知られるNURBSは、Non Uniform Rational B-Splinesの略で、高度なSpline補間方式が用いられています。
線形補間
相隣る二点間を直線（一次式）で結ぶ一番単純な補間で、学校のレポートでグラフを描く時には、「最もいけない」とされるものです。プログラム的にも面白味がありませんので、ここでは取り上げません。
Lagrange補間
データ点を（X0,Y0）、（X1,Y1）、…、（Xn,Yn）の（n+1）個とすると、これらの点を必ず通るn次の多項式Pn（x）は下式で求まります。LはLagrange、PはPolynomial（多項式）の意味です。
プログラムでは、データの個数をdata_countとしていますので、jは0からdata_count-1の間で変化させ、
(x-dataX[j])/(dataX[i]-dataX[j])
を掛け合わせて結果をi番目のlagとし、yはlag*dataY[i]を加算して求めます。
Lagrange補間の改良
上記のLagrange補間は、両端付近で大きな異常振動が生じる場合があります。これを緩和するために、次の改良を行います。
- データ点の（X0,Y0）、（X1,Y1）、…、（Xn,Yn）の（n+1）個に対して、両端に合計2個の点を追加します。入力した（X0,Y0）、（X1,Y1）、…、（Xn,Yn）を、（X1,Y1）、（X2,Y2）、…、（Xn,Yn）、（Xn+1,Xn+1）のように名前を付け替え、両端に新たに（X0,Y0）と（Xn+2,Yn+2）を追加します。プログラムでは、名前の付け替えは、
- （X0,Y0）は、（X1,Y1）と（X2,Y2）を線形的に外挿して求め、（Xn+2,Yn+2）は、（Xn,Yn）と（Xn+1,Yn+1）を線形的に外挿します。
- 以上の（n+3）個の点に対して通常のLagrange補間を行いますが、追加した（X0,Y0）と（X1,Y1）の間、および（Xn,Yn）と（Xn+1,Yn+1）の間は表示しません。
dataX1[i+1]=dataX[i] と dataY1[i+1]=dataY[i]
dataX1[0]=2*dataX[0]-dataX[1] dataX1[data_count+1]=2*dataX[data_count-1]-dataX[data_count-2]
if(x>=dataX1[1] && x<=dataX1[data_count])
実行結果を見ると、ある程度の効果があることが分りますが、線形的な外挿が必ずしも好ましくない場合があり、かえって振動を増加させることがあります。
Newton補間
データ点を（X0,Y0）、（X1,Y1）、…、（Xn,Yn）の（n+1）個とすると、これらの点を必ず通るn次の多項式Pn（x）は下式で得られます。Ciは差分商（Divided Difference） と呼ばれ、〔XiXi-1....X0〕で表わしますが、この計算は、再帰プログラムで容易に実現できます。
プログラムでは、n[i]は（x-dataX[j]）の乗算（ただし、jは0からi-1の間）で求めます。
c[i]は、メソッドcalculateDividedDiff(i,0)で求め、yは、c[i]*n[i]の加算で求めます。
Spline補間
Spline補間は、末端部の条件や区間の取り方により何種類かありますが、結果はほとんど同じになります。データ点（Xi-1,Yi-1）と（Xi,Yi）の間を結ぶSpline補間の一例であり、広く用いられている3次Splineの式を下記に示します。
ただし、
プログラムでは、3次Splineの第1項から第4項までを、yy0、yy1、yy2、yy3として計算し、加算します。
一次微分はdif1[i]、二次微分はdif2[i]で表わし、上式をそのままコード化します。
プログラム
import java.applet.Applet; import java.awt.*; import java.awt.event.*; public class Interpolation extends Applet implements MouseListener,MouseMotionListener,ActionListener{ static final int X0=50,Y0=50; //「方眼紙」の左上座標 int data_count; //データ個数 //入力データのｘ値（座標値ではない） double[] dataX=new double[20]; //入力データのｙ値（座標値ではない） double[] dataY=new double[20]; Button[] bt=new Button[6]; public void init(){ data_count=0; // ---------------------- マウス関係準備 ---------------------- addMouseMotionListener(this); addMouseListener(this); // ---------------------- ボタン関係準備 ---------------------- String[] label={" Lagrange "," Newton ", " Lagrange改 "," ３次Spline ", " 補間消去 "," データ消去 "}; for(int i=0;i<6;i++){ bt[i]=new Button(label[i]); add(bt[i]); bt[i].addActionListener(this); } //ボタンの色を設定する bt[0].setBackground(Color.green); //Lagrange bt[1].setBackground(Color.yellow); //Newton bt[2].setBackground(Color.magenta); //Lagrange改 bt[3].setBackground(Color.cyan); //３次Spline } // ------------------- マウス移動関係メソッド -------------------- public void mouseMoved(MouseEvent me){ //マウスを動かした int x=me.getX(); int y=me.getY(); float xx=(float)((x-X0)/10.0); //floatにキャストして丸める float yy=(float)(40.0-(y-Y0)/10.0); Graphics g=getGraphics(); g.setColor(Color.white); g.fillRect(X0-5,Y0+445,300,20); //前の表示を消す g.setColor(Color.black); if(xx>=0.0 && xx<=60.0 && yy>=0.0 && yy<=40.0) // 座標の表示 g.drawString("X座標= "+ xx +" ,Y座標= "+ yy ,X0,Y0+460); else g.drawString("範囲外です。",X0,Y0+460); } //使用しないが、記述しておかないとエラーになる public void mouseDragged(MouseEvent me){ } // ------------------- マウス動作関係メソッド -------------------- public void mousePressed(MouseEvent me){ //マウスをクリックした int x=me.getX(); int y=me.getY(); double xx=(x-X0)/10.0; double yy=40.0-(y-Y0)/10.0; if(xx>=0.0 && xx<=60.0 && yy>=0.0 && yy<=40.0){ dataX[data_count]=xx; dataY[data_count]=yy; data_count++; repaint(); } } //使用しないが、記述しておかないとエラーになる public void mouseClicked(MouseEvent me){ } public void mouseEntered(MouseEvent me){ } public void mouseReleased(MouseEvent me){ } public void mouseExited(MouseEvent me){ } // --------------------- ボタン関係メソッド ---------------------- //ボタンがクリックされた public void actionPerformed(ActionEvent ae){ if(ae.getSource()==bt[0]) drawLagrange(); //Lagrange if(ae.getSource()==bt[1]) drawNewton(); //Newton if(ae.getSource()==bt[2]) drawLagrange1(); //Lagrange改 if(ae.getSource()==bt[3]) drawSpline(); //３次Spline if(ae.getSource()==bt[4]) repaint(); //補間消去 if(ae.getSource()==bt[5]) clearData(); //データ消去 } // --------------------- 補間処理メソッド類 ---------------------- //Lagrange補間 private void drawLagrange(){ double x,y,lag; Graphics g=getGraphics(); g.setColor(new Color(0,192,0)); //dark greenで線を引く for(x=dataX[0];x<=dataX[data_count-1];x+=0.01){ y=0.0; for(int i=0;i<data_count;i++){ lag=1.0; for(int j=0;j<data_count;j++) if(j!=i) lag*=(x-dataX[j])/(dataX[i]-dataX[j]); y+=lag*dataY[i]; } g.drawRect(X0+(int)(10*x),Y0+(int)(10*(40.0-y)),0,0); } } //Newton補間 private void drawNewton(){ int i,j; double x,y,p; double[] n=new double[20]; double[] c=new double[20]; Graphics g=getGraphics(); g.setColor(new Color(192,192,0)); //dark yellowで線を引く for(x=dataX[0];x<=dataX[data_count-1];x+=0.01){ n[0]=1.0; for(i=1;i<data_count;i++){ p=1.0; for(j=0;j<i;j++) p*=x-dataX[j]; n[i]=p; } for(i=0;i<data_count;i++) c[i]=calculateDividedDiff(i,0); y=0.0; for(i=0;i<data_count;i++) y+=c[i]*n[i]; g.drawRect(X0+(int)(10*x),Y0+(int)(10*(40.0-y)),0,0); } } //Newton補間から呼び出される再帰メソッド private double calculateDividedDiff(int a,int b){ if(a<=b) return dataY[a]; else return (calculateDividedDiff(a,b+1)- calculateDividedDiff(a-1,b))/(dataX[a]-dataX[b]); } //Lagrange補間の改良 private void drawLagrange1(){ double x,y,lag; double[] dataX1=new double[22]; //データ点を増やす準備 double[] dataY1=new double[22]; for(int i=0;i<data_count;i++){ //入力したデータ点を移す dataX1[i+1]=dataX[i]; dataY1[i+1]=dataY[i]; } //両端のデータを外挿して求める dataX1[0]=2*dataX[0]-dataX[1]; dataY1[0]=2*dataY[0]-dataY[1]; dataX1[data_count+1]=2*dataX[data_count-1]-dataX[data_count-2]; dataY1[data_count+1]=2*dataY[data_count-1]-dataY[data_count-2]; Graphics g=getGraphics(); g.setColor(new Color(192,0,192)); //dark magentaで線を引く for(x=dataX1[0];x<=dataX1[data_count+1];x+=0.01){ y=0.0; for(int i=0;i<data_count+1;i++){ lag=1.0; for(int j=0;j<data_count+1;j++) if(j!=i) lag*=((x-dataX1[j])/(dataX1[i]-dataX1[j])); y+=(dataY1[i]*lag); } //実際に入力した範囲のみを描画 if(x>=dataX1[1] && x<=dataX1[data_count]) g.drawRect(X0+(int)(10*x),Y0+(int)(10*(40.0-y)),0,0); } } //３次Spline補間 private void drawSpline(){ int i; double x,y,yy0,yy1,yy2,yy3; double[] h=new double[20]; //間隔 double[] dif1=new double[20]; //一次微分 double[] dif2=new double[20]; //二次微分 h[0]=0.0; dif2[0]=0.0; dif2[data_count-1]=0.0; Graphics g=getGraphics(); g.setColor(new Color(0,192,192)); //dark cyanで線を引く for(i=1;i<data_count;i++){ h[i]=dataX[i]-dataX[i-1]; //間隔を計算 dif1[i]=(dataY[i]-dataY[i-1])/h[i]; //一次微分を計算 } for(i=1;i<data_count;i++) //二次微分を計算 dif2[i]=(dif1[i+1]-dif1[i])/(dataX[i+1]-dataX[i-1]); i=1; for(x=dataX[0];x<dataX[data_count-1];x+=0.01){ if(x<dataX[i]){ yy0=dif2[i-1]/(6*h[i])*(dataX[i]-x) *(dataX[i]-x)*(dataX[i]-x); //第１項 yy1=dif2[i]/(6*h[i])*(x-dataX[i-1]) *(x-dataX[i-1])*(x-dataX[i-1]); //第２項 yy2=(dataY[i-1]/h[i]-h[i]*dif2[i-1]/6) *(dataX[i]-x); //第３項 yy3=(dataY[i]/h[i]-h[i]*dif2[i]/6)* (x-dataX[i-1]); //第４項 y=yy0+yy1+yy2+yy3; g.drawRect(X0+(int)(10*x),Y0+(int)(10*(40-y)),0,0); } else i++; } } //データ点消去 private void clearData(){ data_count=0; for(int i=0;i<20;i++){ dataX[i]=0.0; dataY[i]=0.0; } repaint(); } public void paint(Graphics g){ int i,j,x,y; //「方眼紙」を描く g.setColor(Color.lightGray); for(i=0;i<=60;i++) g.drawLine(X0+10*i,Y0,X0+10*i,Y0+400); for(j=0;j<=40;j++) g.drawLine(X0,Y0+10*j,X0+600,Y0+10*j); g.setColor(Color.gray); for(i=0;i<=12;i++) g.drawLine(X0+50*i,Y0,X0+50*i,Y0+400); for(j=0;j<=8;j++) g.drawLine(X0,Y0+50*j,X0+600,Y0+50*j); //目盛を入れる g.setColor(Color.black); for(i=0;i<=12;i++) g.drawString(String.valueOf(5*i),X0-5+50*i,Y0+424); for(j=0;j<=8;j++) g.drawString(String.valueOf(5*(8-j)),X0-25,Y0+5+50*j); //データ点を描く g.setColor(Color.red); i=data_count-1; while(i>=0){ x=X0+(int)(dataX[i]*10); y=Y0+(int)((40.0-dataY[i])*10); g.drawOval(x-2, y-2, 4, 4); i--; } } }
プログラムの使い方
アプレットを実行すると、上部に、左から［Lagrange］［Newton］［Lagrange改］［3次Spline］［補間消去］［データ消去］のボタンが表示され、下に方眼紙状のデータ表示領域が現われます。
「方眼紙」内でクリックすると、データ点が入力され、赤丸が表示されます。マウスポインタの位置は、左下に数字で表わされますので、参考にします。データ入力は、必ず左から行ってください。途中で引き返すと、その点までしか表示されません（これは改善できますが、プログラムを簡単化するために省略しました）。間違ってデータを入力した場合は、［データ消去］のボタンをクリックして、最初からやり直してください。
また、入力できるデータ点の数は20個に制限してありますので、注意してください。
［Lagrange］などのボタンをクリックすると、補間が実行され、それぞれ、ボタンと同じ色で曲線が描かれます。次々とボタンをクリックすると、対応する補間曲線が追加されます。［補間消去］ボタンをクリックすると、データ点を残して補間曲線のみが消えるので、何回でも繰り返すことができます。
データを入力し直すときは、［データ消去］をクリックしてください。
なお、［Lagrange］と［Newton］は全く同じ結果になります。両者を実行できるようにしたのは、その事実の確認のためです。
プログラムの実行結果
下図において、赤丸は、「暴れ」が出やすいように、やや意地悪に入力したデータ点（普通の実験データの場合は、きれいな曲線になります）、緑の線はLagrange、マゼンタはLagrange改良型、シアンは3次Splineによる補間結果です。Newton補間はLagrange補間と結果が同じなので、表示していません。Lagrangeに比べ、Lagrange改良型は暴れ方が減っていることがわかります。3次Splineは、「暴れ」はありませんが、線がスムーズでない部分があります。
「方眼紙」の下は、マウスポインタが指している座標で、この例では、一番右のデータ点を指しています。
まとめ
ここで紹介した補間プログラムは、プログラム的には比較的単純なものです。各種の補間方式の数学的意味が理解できれば、後は、それをコーディングするだけで、特別なアルゴリズムも要りません。プログラムが得意な人は、数学も得意でしょうから（筆者は少なくともそれを期待しています）、数学的な詳しい説明は省きました。
このプログラムは、プログラムを作成する楽しみもありますが、実際にデータを入れてみて、どのような結果になるかを試してみる面白さがあります。電気に詳しい人であれば、補間結果の「暴れ」の原因がお分りかと思います。データ群は高い周波数成分を有するのに、強制的にn次式でしか表わされず、高調波成分がカットされるからです。
参考資料
この記事は、筆者のウエブサイト『Visual C++の易しい使い方－Lagrange、Newton、Splineなどの色々な補間方式を実現する－（ただし、現在はVisual C++ 2005 Express Edition版に改定）』をJava言語に書き直し、分かりやすく加筆したものです。
補間は、数値計算の基本で、初期のコンピュータの技術計算への応用の主なものの一つでした。したがって、
- 『岩波数学辞典第3版』 日本数学会 編、岩波書店、1985年12月
などの数学辞典類にはたいてい載っています。プログラムに限れば、
- 『C & FORTRANによる数値解析の基礎』 川崎晴久 著、共立出版、1993年5月
- 『CとJavaで学ぶ数値シミュレーション入門』 峯村吉泰 著、森北出版、1999年4月
- 『ANSI Cによる数値計算法入門』 堀之内総一 他著、森北出版、2002年12月
などがあります。
