2010.10.30.
2014. 2.21.訂正 double → float
石立 喬
Visual C++ 2010 Express の易しい使い方（１7）
―――― Lagrange、Newton、Splineなどの色々な補間方式を実現する ――――
概要
実験結果をグラフに表示する時、測定点をプロットしてから、その間を自在定規で結んだことがあると思う。これを自動的に行うプログラムを紹介する。この作業は「補間」と呼ばれる。補間の方式は数多く知られ、プログラムは比較的簡単なので、ボタンを付けてWindowsプログラムらしくした。
補間
「狭義の補間（Interpolation）」は、二つの点におけるデータを絶対に正確なものと仮定し、その中間点のデータを推測することを言う。点が実験などの測定点である場合には、測定データに誤差があっても、そのまま補間を行うので、補間はより正しい値を推測するものではない。「広義の補間」は、与えられたすべての点を通るように曲線を描くこと、その曲線の方程式を求めることである。
補間のいろいろ
一番単純な補間は二点間を直線（一次式）で結ぶもので、「線形補間」と呼ばれる。学校のレポートでグラフを描く時には、「最もいけない」とされるものである。0番目からn番目まで（n+1）個の点を必ず通る曲線（測定点のデータに誤差があっても必ず通るのは良くないが）を表現するには、ｎ次の多項式が必要十分条件であるので、この多項式を求めることを一般に補間と呼ぶ。
多項式を求める方法には、Lagrange多項式が良く知られているが、再帰プログラムを用いるNewtonの方法もあり、二つの方法による結果は同じになる。ここで紹介するLagrange補間の改良は、筆者の案によるものである（すでに誰かがやっているかも知れないが…）。一方、多項式の次数を限定し、狭い区間毎に異なる多項式を求める方法にSpline法がある。
Lagrange補間
データ点を(X0,X0)、(X1,Y1)、…、(Xn,Yn)の（n+1）個とすると、これらの点を必ず通るｎ次の多項式Pn(x)は下式で求まる。LはLagrange、PはPolynomial（多項式）の意。
Lagrange補間の改良
上記のLagrange補間は、両端で大きな振動が生じる場合がある。これを緩和するために、データ点の(X0,Y0)、(X1,Y1)、…、(Xn,Yn)の（n+1）個に対して、両端に２個の点を追加する。入力した(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)の間は表示しない。データ点の条件にもよるが、実行結果を見ると、ある程度の効果があることが分る。
Newton補間
データ点を(X0,Y0)、(X1,Y1)、…、(Xn,Yn)の（n+1）個とすると、これらの点を必ず通るｎ次の多項式Pn(x)は下式で求まる。〔XiXi-1
… X0〕の計算は、再帰プログラムで容易に実現できる。
Spline補間
Spline補間は、末端部の条件や区間の取り方により何種類かあるが、結果はほとんど同じである。データ点(Xi-1,Yi-1)と(Xi,Yi)の間を結ぶSpline補間の式の一例（最も良く使用される3次Spline補間）を下記に示す。
ただし、
プログラムの概要
◎準備
１）マウスで入力するデータ点は、多くのメソッドから呼び出されるので、 配列xdata[i]、ydata[i]をそれぞれ、20個だけグローバルに設定する。staticにする必要がある。
２）Form1_Paint()では、グラフ用紙を作成し、目盛を入れ、データ点があれば、それらを赤丸で表示する。
◎マウス操作
１）マウスをグラフ用紙の領域内で移動させると、その時のグラフ用紙上の値を、画面左下に表示する。グラフ用紙上の値範囲外の場合には、「Out
of Range」と表示する。これは、Form1_MouseMove()を使用する。マウスによって指定されたパソコン画面上の絶対位置座標（e->X,e->Y）は、
x=(e->X-X0)/10;
y=40-(e->Y-Y0)/10;
により、グラフ用紙上の値に変換される。ただし、xおよびyは、float型とする。
２）マウスをグラフ用紙の領域内でクリックすると、その時のグラフ用紙上の値を、配列xdata[i]、ydata[i]に格納し、data_numberをインクリメントし、Invalidate()でForm1_Paintを呼び出して、クリックした場所を赤丸で表示する。最初の値はx_data[0],y_data[0]に、最後の値はx_data[data_number-1],y_data[data_number-1]に入る。
◎各補間プログラム
１）共通事項
グラフ用紙上で、ｘをxdata[0]からxdata[data_number-1]まで、0.2f刻みで変化させ、各補間式を実行し、式によって得たｘとｙの値を、
xx=X0+x*10
yy=Y0+(40-y)*10
によってパソコン画面上の絶対位置座標に変換する。
グラフは、点（xx_old,yy_old）から点（xx,yy）までDrawLineによって描画する。「Lagrange改」を除いて、x=xdata[0]の時は、xx_oldとyy_oldを設定するだけで、線は引かない。「Lagrange改」の場合には、指定範囲外では線を引かないので、flagを用いて制御する。
２）Lagrange補間（buttonLagrange_Click）
Lagrange補間の式をそのまま忠実に実行する。乗算を繰り返して行なうので最初にlag=1.0f;
と初期設定しておく。
３）Newton補間（buttonNewton_Click）
配列nとcを、それぞれ20個用意する。これは、入力可能なデータ個数の上限に合わせてある。c[i]を求めるために、自作のメソッドdifference(a,b)を用意する。このメソッドは、前述のNewton補間の式に従って作られていて、再帰的に自分を呼び出す。
４）Lagrange補間の改良（buttonLagrange1_Click）
データ点に対して両端の2点を付加するので、新しくデータ点xdata1とydata1の配列をそれぞれ22個用意する。
実際に入力したデータ点をxdata1[i+1]=xdata[i];
などにより一つずつずらし、両端、すなわち、xdata1[0]とxtata1[data\number+1]などを線形的に外挿して求める。
新しく作ったデータ点について通常のLagrange補間の計算をした後、実際に入力したデータ点の間だけに線を描く。初めて線を描く場合にはold_xxやold_yyを設定するのみにとどめるため、flagを用いて制御する。
５）３次Spline補間（buttonSpline_Click）
Spline補間の式にしたがって、配列h（データ間隔）、dif1（一次微係数）、dif2（二次微係数）をそれぞれ20個用意し、h[0]などを0.0fに設定する。後は、式にしたがって計算する。
◎ボタン関係
１）ボタン[補間消去]をクリックした場合(buttonInterCancel_Click)
Invalidate()でForm1_Paintを呼び出して、全てのグラフィック描画を一旦消去し、データ点のみを赤丸で表示する。補間結果は描画しないので、補間結果は消去される。
２）ボタン[データ消去]をクリックした場合(buttonDataCancel_Click)
データの個数data_numberをゼロにクリヤし、配列xdata[i]、ydata[i]もゼロに設定する。次いで、Invalidate()でForm1_Paintを呼び出して、全てのグラフィック描画を消去する。データ点の個数がゼロなので、赤丸は表示されない。
実数にはdoubleではなくfloatを使う
一般的には、整数はint、実数はdoubleを使うのが原則である。コンピュータの高速化、メモリの大容量化により、shortやfloatを使う必要はなくなった。
ところが、DrawLineメソッドやDrawEllipseメソッドは、座標をint（Int32）またはfloat（Single）で与えることになっていて（他にRectangleやRectangleFも可）、doubleで与えると警告が、doubleとintの混在で与えるとエラーが表示される。
実数をdoubleで計算し、DrawLineなどのメソッドにはfloatにキャストして使用するのは面倒なので、ここでは最初からfloatで計算した。
プログラム
static int X0=50,Y0=60;
int data_number;
static array<float>^ xdata=gcnew array<float>(20);
static array<float>^ ydata=gcnew array<float>(20);
private: System::Void Form1_Paint(System::Object^ sender, System::Windows::Forms::PaintEventArgs^ e) {
int i,j;
float xx,yy;
Graphics^ g=e->Graphics;
//薄く線を引く
for(i=0;i<=60;i++)
g->DrawLine(Pens::LightGray,X0+10*i,Y0,X0+10*i,Y0+400);
for(j=0;j<=40;j++)
g->DrawLine(Pens::LightGray,X0,Y0+10*j,X0+600,Y0+10*j);
//濃く線を引く
for(i=0;i<=12;i++)
g->DrawLine(Pens::Gray,X0+50*i,Y0,X0+50*i,Y0+400);
for(j=0;j<=8;j++)
g->DrawLine(Pens::Gray,X0,Y0+50*j,X0+600,Y0+50*j);
//目盛を入れる
String^ string1;
for(i=0;i<=12;i++){
string1=String::Format(" {0} ",5*i);
g->DrawString(string1,Font,Brushes::Black,X0-15+50*i,Y0+410);
}
for(j=0;j<=8;j++){
string1=String::Format(" {0} ",5*(8-j));
g->DrawString(string1,Font,Brushes::Black,X0-30,Y0-6+50*j);
}
//データ点を描画する
for(i=0;i<data_number;i++){
xx=X0+xdata[i]*10;
yy=Y0+(40-ydata[i])*10;
g->DrawEllipse(Pens::Red,xx-2,yy-2,4.0f,4.0f);
}
}
private: System::Void Form1_MouseMove(System::Object^ sender, System::Windows::Forms::MouseEventArgs^
e) {
Graphics^ g=this->CreateGraphics();
//データ値に変換する
float x=(e->X-X0)/10.0f;
float y=40-(e->Y-Y0)/10.0f;
String^ string1;
//前の表示を消して、新しいデータ値を表示する
g->FillRectangle(Brushes::White,X0-5,Y0+446,90,20);
if(x>=0.0f && x<=60.0f && y>=0.0f && y<=40.0f)
string1=String::Format("X={0}, Y={1}",x,y);
else
string1="Out of Range";
g->DrawString(string1,Font,Brushes::Black,X0,Y0+450);
}
private: System::Void Form1_MouseDown(System::Object^ sender, System::Windows::Forms::MouseEventArgs^
e) {
Graphics^ g=this->CreateGraphics();
//データ値に変換する
float x=(e->X-X0)/10.0f;
float y=40-(e->Y-Y0)/10.0f;
//データ値を格納する
if(x>=0.0f && x<=60.0f && y>=0.0f && y<=40.0f){
xdata[data_number]=x;
ydata[data_number]=y;
data_number++;
}
Invalidate();
}
private: System::Void buttonInterCancel_Click(System::Object^ sender, System::EventArgs^ e) {
Invalidate();
}
private: System::Void buttonDataCancel_Click(System::Object^ sender, System::EventArgs^
e) {
//データ個数をリセットする
data_number=0;
//データ値をリセットする
for(int i=0;i<20;i++){
xdata[i]=0.0f;
ydata[i]=0.0f;
}
Invalidate();
}
private: System::Void buttonLaglange_Click(System::Object^ sender, System::EventArgs^
e) {
Graphics^ g=this->CreateGraphics();
int i,j;
float x,y,lag;
float xx,yy,xx_old,yy_old;
for(x=xdata[0];x<=xdata[data_number-1];x+=0.2f){
y=0.0f;
for(i=0;i<data_number;i++){
lag=1.0f;
for(j=0;j<data_number;j++)
if(j!=i) lag*=(x-xdata[j])/(xdata[i]-xdata[j]);
y+=ydata[i]*lag;
}
//画面上の座標に変換する
xx=X0+x*10;
yy=Y0+(40-y)*10;
if(x==xdata[0]){
xx_old=xx;
yy_old=yy;
}
else{
g->DrawLine(Pens::Green,xx_old,yy_old,xx,yy);
xx_old=xx;
yy_old=yy;
}
}
}
private: System::Void buttonNewton_Click(System::Object^ sender, System::EventArgs^
e) {
Graphics^ g=this->CreateGraphics();
int i,j;
float x,y,p;
float xx,yy,xx_old,yy_old;
array<float>^ n=gcnew array<float>(20);
array<float>^ c=gcnew array<float>(20);
for(x=xdata[0];x<=xdata[data_number-1];x+=0.2){
n[0]=1.0f;
for(i=1;i<data_number;i++){
p=1.0f;
for(j=0;j<i;j++)
p*=x-xdata[j];
n[i]=p;
}
for(i=0;i<data_number;i++)
c[i]=difference(i,0);
y=0.0f;
for(i=0;i<data_number;i++)
y+=c[i]*n[i];
xx=X0+x*10;
yy=Y0+(40-y)*10;
if(x==xdata[0]){
xx_old=xx;
yy_old=yy;
}
else{
g->DrawLine(Pens::DarkOrange,xx_old,yy_old,xx,yy);
xx_old=xx;
yy_old=yy;
}
}
}
private: double difference(int a,int b) {
if(a<=b) return ydata[a];
else return (difference(a,b+1)-difference(a-1,b))/(xdata[a]-xdata[b]);
}
private: System::Void buttonLaglange1_Click(System::Object^ sender, System::EventArgs^ e) {
Graphics^ g=this->CreateGraphics();
int i,j;
float x,y,lag;
float xx,yy,xx_old,yy_old;
array<float>^ xdata1=gcnew array<float>(22);
array<float>^ ydata1=gcnew array<float>(22);
bool flag=true;
//入力したデータ点を一つずらす
for(i=0;i<data_number;i++){
xdata1[i+1]=xdata[i];
ydata1[i+1]=ydata[i];
}
//入力したデータ点を外挿する
xdata1[0]=2*xdata[0]-xdata[1];
ydata1[0]=2*ydata[0]-ydata[1];
xdata1[data_number+1]=2*xdata[data_number-1]-xdata[data_number-2];
ydata1[data_number+1]=2*ydata[data_number-1]-ydata[data_number-2];
for(x=xdata1[0];x<=xdata1[data_number+1];x+=0.2f){ //外挿データ点を含む
y=0.0f;
for(i=0;i<data_number+2;i++){
lag=1.0f;
for(j=0;j<data_number+2;j++)
if(j!=i) lag*=(x-xdata1[j])/(xdata1[i]-xdata1[j]);
y+=ydata1[i]*lag;
}
if(x>=xdata[0] && x<=xdata[data_number-1]){ //実データ点のみ描画
xx=X0+x*10;
yy=Y0+(40-y)*10;
if(flag){ //描画する最初であったら
xx_old=xx;
yy_old=yy;
flag=false;
}
else{
g->DrawLine(Pens::Magenta,xx_old,yy_old,xx,yy);
xx_old=xx;
yy_old=yy;
}
}
}
}
private: System::Void buttonSpline_Click(System::Object^ sender, System::EventArgs^
e) {
Graphics^ g=this->CreateGraphics();
int i;
float x,y,yy0,yy1,yy2,yy3;
float xx,yy,xx_old,yy_old;
array<float>^ h=gcnew array<float>(20);
array<float>^ dif1=gcnew array<float>(20);
array<float>^ dif2=gcnew array<float>(20);
h[0]=0.0f;
dif2[0]=0.0f;
dif2[data_number-1]=0.0f;
for(i=1;i<data_number;i++){
h[i]=xdata[i]-xdata[i-1];
dif1[i]=(ydata[i]-ydata[i-1])/(xdata[i]-xdata[i-1]);
}
for(i=1;i<data_number-1;i++)
dif2[i]=(dif1[i+1]-dif1[i])/(xdata[i+1]-xdata[i-1]);
i=1;
for(x=xdata[0];x<xdata[data_number-1];x+=0.2){
if(x<xdata[i]){
yy0=dif2[i-1]/(6*h[i])*(xdata[i]-x)*(xdata[i]-x)*(xdata[i]-x);
yy1=dif2[i]/(6*h[i])*(x-xdata[i-1])*(x-xdata[i-1])*(x-xdata[i-1]);
yy2=(ydata[i-1]/h[i]-h[i]*dif2[i-1]/6)*(xdata[i]-x);
yy3=(ydata[i]/h[i]-h[i]*dif2[i]/6)*(x-xdata[i-1]);
y=yy0+yy1+yy2+yy3;
xx=X0+x*10;
yy=Y0+(40-y)*10;
if(x==xdata[0]){
xx_old=xx;
yy_old=yy;
}
else{
g->DrawLine(Pens::DarkCyan,xx_old,yy_old,xx,yy);
xx_old=xx;
yy_old=yy;
}
}
else i++;
}
}
補間プログラムの実行結果の一例
下図において、赤丸は適当にマウスで入力したデータ点、緑の線はLagrange補間、マゼンタはLagrange補間の改良型、シアンは３次Spline補間によるものである。Newton補間の結果はLagrange補間の結果と全く同じなので、ここでは表示していない。Lagrangeに比べ、Lagrange改良型は両端における暴れ方が減っていることが分る。３次Splineは、暴れはないが、線がスムーズでない部分がある。「方眼紙」の下は、マウスが指している座標である（この例では、一番右のデータ点を指している）。 このプログラムは、プログラムを作成する楽しみもあるが、実際にデータを入れてみて、どのような結果になるかを試してみる面白さもある。電気に詳しい諸君であれば、補間方程式で得られた結果の「暴れ」の原因がお分かりと思う。データ群は高い周波数成分を有するのに、強制的にｎ次式でしか表されず、高調波成分がカットされるからである。
付録：ボタンの作り方
「易しい使い方（３）」で、すでにボタンの作り方について簡単に述べてあるが、あらためて解説する。
１）「Form1.h[デザイン]」を画面に表示する。
２）「表示」→「ツールボックス」をクリックする。
３）ツールボックスで、「Button」コントロールをクリックする。
４）四角い枠に「ab」と書かれたカーソルが、マウスに付いてくるので、Form1上の希望する位置でクリックする。大体の位置やサイズを決めておいて、詳細は、「Form1.h」の上部にある
this->buttonLaglange->Location = System::Drawing::Point(40,
18);
this->buttonLaglange->Size = System::Drawing::Size(75, 23);
などの括弧内を直接編集して揃えるとよい（「コード エディターで変更しないでください。」とあるが、大丈夫）。
５）設定したいボタンを右クリックし、「プロパティ」ウインドウで、
デザイン ---- (Name) 「button1」などを、「buttonLanguage」など分かりやすい名前に変更
表示 ------- BackColor 右にある「▼」をクリックし、「カスタム」、「Web」、「システム」タブを切り替えて選択
「LightGreen」などと直接入力も可
Text 「button1」などを、「Laglange」などに変更
６）各ボタンをダブルクリックすると、イベントハンドラメソッドのスケルトンが現れるので、そこに内容を記述する。
「buttonClear」などの名称を間違えて設定してしまった場合には、
１）ボタンを右クリックして、「削除」をクリックする。
２）「Form1.h」で、関連するイベントハンドラ関数を削除する。
の手続きが必要である。
「Visual C++ の勉強部屋」（目次）へ