第2回　補間とスプライン(1)

本格的な工学計算に向けて，要点を確認しておこう．

「配列の構造体と，「構造体の配列」

空間座標やベクトルなどの量を扱う場合は，複数の値をひとまとめにして考えることが多い．例えば，ベクトルでは二つや三つの実数を組み合わせて1つのベクトルが定義されるし，二次元平面内の任意の点は二つの実数 (x, y)で表される．
また，複素数は実部と虚部の2つの実数の組み合わせで表される．（C言語では複素数を直接的に扱う組み込み変数は用意されていないため，自分で用意するか，または複素数用の外部の拡張ライブラリを使用する必要がある．）
C言語でこのような概念を表現する場合，構造体が用いられる．

以下に，例として10個の二次元座標 ( x, y ) の表現方法を考える．

( x, y ) の配列をそれぞれ用意する方法．


int main()
{
    double x[10];
    double y[10];

    x[0] = 1.0;   /*  同じ添字が　単一の点を表す．  */
    y[0] = 2.0;

    ...
}

配列を構造体のメンバ変数とする方法．


struct Points2Darray
{
    double x[10];
    double y[10];
};

int main()
{
    Point2Darray points;

    points.x[0] = 1.0;   /*  構造体の中身（メンバ変数）が配列  */
    points.y[0] = 2.0;
    ...
}

構造体を配列とする方法．


struct Point2D
{
    double x;
    double y;
};

int main()
{
    Point2D points[10];

    points[0].x = 1.0;    /*  配列の中身（各要素）が構造体  */
    points[0].y = 2.0;
	...
}

いずれの方法も，10個の二次元座標(x,y)を正しく取り扱うことができるので，どれを用いても構わないが，

配列へのアクセス方法が異なる．（上記の main 関数内を参照）
メモリ上のデータの並び方が，1, 2の場合はx[0], x[1], ... , x[99], y[0], y[1], ... , y[99]となるが，3ではx[0], y[0], x[1], y[1], ... , x[99], y[99]となる．
前者は SoA (Structure of Array)，後者は AoS (Array of Structure) と呼ばれる．
関数に構造体を渡すと，値参照（変数全部がコピーされる）なので非効率．1,3 は配列なのでアドレス参照（先頭アドレスのみ）ので無駄がない．もちろん前者でもアドレスで渡せば効率的．
これらは処理しようとする計算に応じて向き不向きがあり，その計算する目的に応じて使い分けられる．一般に，一度の計算に使用する変数同士がメモリ上の近くに有る方が処理速度が速くなる．
特にSuperComputer ( =並列計算機) を用いた大規模計算においては，この差は極めて重要となる．

補間

コンピュータは有限個のデータ（離散データという）であれば処理することができるが，値が無限に必要な連続関数を数値データとして扱うことはできない．
そこで，有限個の離散データをもとに，まずその点を必ず通る関数の形（具体的には係数）を決定し，それを用いてデータとして持っていない任意の関数値を得る方法を補間と呼ぶ．（「補完」ではないので注意．英語では interpolation）

（これに対し，与えた点群を通らなくても良く，ある関数に対して点群との誤差が最小となるように係数を決定する手法を最小二乗近似と呼ぶ．後半の実習で取り上げる．）

補間には多くの種類があるが，幾つか基本的な方法を記す．

線形補間 (linear interpolation) ，直線補間とも言う．
ラグランジュ補間 (Lagrange interpolation)
スプライン補間 (spline interpolation)

図1. 線形補間の例．既知の2点 (x₀, y₀), (x₁, y₁) から1次関数が決定できれば，任意のxにたいする関数値yを求めることができる．

もっとも単純な補間法は線形補間である．線形補間では，2点が与えられればそれを直線的に結ぶ関数を決定することができるため，その区間内外の任意の点に対する関数値を決定することができる．

線形補間の式

任意の2点に対する線形補間の式は以下のように表される．

y - y₀ = (( y₁ - y₀) ⁄ (x₁ - x₀)) (x - x₀)　　　　　但し， x₀≠x₁ 　　　　　　　 (1)

これを用いて，任意のxにたいするyを計算することができる．
より高次の関数についても，多項式の各係数を決めることにより関数が決定できる．

内挿と外挿

点群から補間により関数を決定し，任意の場所の関数値を計算できるが，どの位置の値を得るかにより内挿と外挿とに分けられる．

図2. 内挿と外挿の比較．点群の両端の区間内の点（ x₀ ≤ x ≤ x_N ）の点を求めることを内挿, interpolation，区間外の点を計算することを外挿, extrapolationと言う．

内挿 (interpolation)

午前9時と10時の気温から，9時30分の気温を推測．
幾つかの制御点のみでなめらかな曲線を表現．（ゲームなどのCG，画面やプリンタへの文字表示）
音声データのレート変換処理．
画像の拡大・縮小，回転などの処理．

外挿 (extrapolation)

今日までの株価を元に，明日の株価を予想．
今日までの過去の気温を元に，明日の気温を予想．

一般的に内挿に比べて外挿はその値の信頼性について担保することが難しいことが知られている．

課題1

以下の手順に従ってソースファイルおよびグラフ等を作成せよ．

(1)

点群の座標値を画面に表示する関数disp_xy()を作成せよ．テスト用main関数および実行結果を以下の通りとせよ．


struct Point2D
{
    double x;
    double y;
};

void disp_xy(const Point2D* pnt, const int n)  /* この関数内で仮引数を変更しない場合はこのように書く．*/
{
    ...
}

int main()
{
    const int N = 5;
    Point2D points[N] =
       {{-3.5, -2.2},
        {-2.2, -0.3},
        {-1.0,  1.5},
        { 0.2,  1.9},
        { 2.3,  2.2}};

    disp_xy(points, N);
}

実行結果：
-3.500000 , -2.200000
-2.200000 , -0.300000
-1.000000 , 1.500000
0.200000 , 1.900000
2.300000 , 2.200000

(2)

disp_xy()を改良して，点群の座標をcsvファイル 153Rxxxxxx-03-1.csv に出力(xxxxxxxは各自の学生番号)する関数fileout()を作成せよ．
引数および戻り値はdisp_xy()と同じで良い．
csvファイル内の各行に，カンマ区切りでx,yの座標値が格納されるようにすること．

(3)

前問で得られた csv ファイルをExcelで開いて.xlsx形式で保存し，散布図としてグラフを書け．
完成した.xlsxファイルおよびソースファイルを提出せよ．

散布図の例

課題2

前問において，各区間（計4区間）について，2点を通る1次関数を，式(1)を用いて求めよ．計算結果として，x⁰およびx¹の係数（つまり「切片」と「傾き」）を画面に表示せよ．ソースファイルを提出すること．

可能であれば，計算結果が正しいかどうか，グラフ用紙上，またはExcelで検証すること．（得られた切片と傾きから直線を書き，直線上に2点が存在するかを確認．）

実行結果（数値は適当であり，正しくない）

        | 切片　　   傾き
--------+------------------
区間1   | 5.234    0.234
区間2   | 2.300    1.113
区間3   | -1.233   3.224
区間4   | 2.344    1.11

予習１：高次関数による補間法ラグランジュ補間（多項式）

上記課題では，2点から1次関数による補間を行ったが，実際には2次関数，3時間数などの高次の多項式による補間を行う場合も多い．一般に，N+1個の座標値 (x₀, y₀), (x₁, y₁), (x₂, y₂), ... , (x_N, y_N) が与えられた場合に，これらの全ての点を通る関数を求める．

N+1 個の値から，N次関数を1つ決定することができる．
例えば，2個の点が決まれば1次関数（直線），3個の点からは2次関数（放物線）が1つ決まる．

一般に，N+1個のデータからは N次式が決定でき，その際に必要な係数はN+1個となる．

y = a₀ + a₁ x + a₂ x² + ... + a_N x^N 　　　　　(1)

係数 a_iをすべて求めれば関数が決まり，関数が決まれば任意の x に対する yの値が得られる．
一般的にこの係数は，N+1元の連立1次方程式を解くことにより求めることができる．
連立方程式の解法には多くの計算時間を要するが，ここでは良い方法がある．
各自で調べてみよう．

予習２：スプライン補間

前述のラグランジュ補間では，多項式で任意の関数を表すものであるが，いくつか問題がある．例えば，周期関数や指数関数，正規分布など，必ずしも多項式での近似がふさわしくない場合がある．（多項式では，xが正負の無限大でその絶対値が無限大になってしまうが，sin, cos などは有限である．）
そこで，与えられた点群すべてを使って関数を決定するのではなく，局所的な数点のみを使って局所の関数（曲線）を決定し，かつ，その曲線が全体として滑らかに接続されるような補間法が広く利用されている．
その代表例がスプライン補間である．次週に向けて，各自で調査・予習しておくこと．

第2回 補間とスプライン(1)