Legendre多項式で関数近似する

最近プログラムを書く機会があったので、そのまとめ。

Legendre多項式とはについての多項式で、
それぞれがすべて正規直交しているので、フーリエ級数のように
データ列を関数近似するのに使われるようです。
（フーリエ級数ほどお目にかからないけど）

具体的には、こんな関数で近似する。

よく見かける基底の係数の公式は連続関数に対するもので、
離散データ（等間隔データも）に関する公式がなぜか見当たらない。

というわけで、基底の離散系列はGram-Schmidtのアルゴリズムを使って自力で求めることになった。

プログラムの手順としては、

1. の離散系列を生成する。
2. Gram-Schmidtのアルゴリズムで系列を修正し、正規直交基底の系列を得る。
3. と入力系列の内積をとって出力係数とする。

これで元の入力系列は、

で近似できることになる。

ソースコードはこんな感じ。割と間に合わせです。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

double w[32][1024];
double in[1024];
double out[1024];

double prod(double a[], double b[], int n)
{
	double ret = 0.0;
	int i;
	for (i = 0; i < n; i++)
		ret += a[i] * b[i];
	return ret;
}

void init(int n, int d)
{
	int d1, d2, i;
	double m;
	
	// generate x^0, x^1
	for (i = 0; i < n; i++)
	{
		w[0][i] = 1.0;
		w[1][i] = 2.0 *  i / (n - 1.0) - 1.0;
	}
	
	// generate x^i
	for (d1 = 2; d1 < d; d1++)
		for (i = 0; i < n; i++)
			w[d1][i] = w[d1 - 1][i] * w[1][i];
	
	// Gram-Schmidt algorithm
	for (d1 = 0; d1 < d; d1++)
	{
		for (d2 = 0; d2 < d1; d2++)
		{
			m = prod(w[d1], w[d2], n);
			for (i = 0; i < n; i++)
				w[d1][i] -= m * w[d2][i];
		}
		
		m = sqrt(prod(w[d1], w[d1], n));
		for (i = 0; i < n; i++)
			w[d1][i] /= m;
	}
}

int main(int argc, char *argv[])
{
	int n, d, i, j;
	FILE *fi = fopen(argv[1], "r");
	FILE *fo = fopen(argv[2], "w");
	n = atoi(argv[3]);
	d = atoi(argv[4]);
	double m;
	
	init(n, d);	
	
	for (i = 0; i < n; i++)
	{
		fscanf(fi, "%lf", &in[i]);
		out[i] = 0.0;
	}
	
	for (i = 0; i < d; i++)
	{
		m = prod(in, w[i], n);
		for (j = 0; j < n; j++)
			out[j] += m * w[i][j];
	}

	for (i = 0; i < n; i++)
		fprintf(fo, "%.8f\n", out[i]);
	
	fclose(fi);
	fclose(fo);
	
	return 0;
}

init()で基底を生成して、残りの処理はmain()で済ませてる。
prod()は系列同士の内積を求める関数です。

テスト入力として、多項式と正弦波を適当に組み合わせてこんな波形を作った。

で、プログラムに食わせた結果がこれ。

Dってのが近似の最高次数（上のプログラムだとatoi(argv[4]) - 1）で、
とりあえず1次近似から17次近似まで4次飛ばしで計算してみた。
次数を上げると、元の波形にどんどん張り付いてゆくのがわかる。