C言語による数値計算を教えて下さい
「2次関数y = ax2 + bx + c の係数を最小2乗法で求めるための正規方程式を導出せよ．その上でキーボードからデータの組(xi, yi)を
5点入力してa,b,cを求めるプログラムを作成せよ(x1～x5,y1～y5)」
という問題なのですが、a,b,cを求める方程式を行列式の形で以下のように出せました。
|Σxi^4 , Σxi^3 , Σxi^2| |a| |Σxi^2*yi|
|Σxi^3 , Σxi^2 , Σxi |*|b| = |Σxi*yi |
|Σxi^2 , Σxi , 5 | |c| |Σyi |
ですがこれをガウス消去法で解こうと思っているのですがどうにもプログラムができません。
どなたかC言語で数値計算できる方、助けてください。
以下のはサイエンス社の「C言語による数値計算入門」をもとに
ガウス消去法で解くプログラムを自分なりに作ったのですが
当然エラーだらけです。
プログラムを作れる方、どうか教えてください。
#include <stdio.h>
#include <math.h>
/* ガウス消去法 */
double *simple_gauss( double a[3][3], double b[3] );
main(void)
{
double a[3][3]; /* マトリックスの宣言*/
double x[3]; /* 解ベクトルの宣言*/
double y[3]; /* 右辺ベクトルの宣言*/
double xx, yy; /* キーボード入力のための一時変数*/
int i, j; /* 繰り返しのための一時変数*/
for(i = 0; i < 2; i++)
{ /* 変数の初期化*/
x[i] = y[i] = 0;
for(j = 0; j < 3; j++)
a[i][j] = 0;
}
for(i=0; i<5; i++)
{
printf("\n(x, y) = ");
scanf("%lf,%lf", &xx, &yy);
a[1][1] += xx * xx * xx * xx;
a[1][3] += xx * xx;
a[2][2] = a[3][1]=a[1][3];
a[1][2] += xx * xx * xx;
a[2][1]=a[1][2];
a[2][3] += xx;
a[3][2]=a[2][3];
a[3][3] ++;
y[2] += xx * yy;
y[1] += xx * xx * yy;
y[3] += yy;
}
y=*simple_gauss(a,y);
return 0;
}
double *simple_gauss( double a[3][3], double y[3] )
{
int i, j, k;
double alpha, tmp;
/* 前進消去 */
for( k = 0; k <= 1; k++)
{
for( i = k+1; i <= ; i++)
{
alpha = - a[i][k]/a[k][k];
for( j = k+1; j <= 2; j++)
{
a[i][j] = a[i][j] + alpha * a[k][j];
}
y[i] = y[i] + alpha * y[k];
}
}
/* 後退代入 */
y[2] = y[2]/a[2][2];
for( k = 1; k >= 0; k--)
{
tmp = y[k];
for( j = k+1; j <= 2; j++)
{
tmp = tmp - a[k][j] * y[j];
}
y[k] = tmp/a[k][k];
}
return y;
}
お礼
ありがとうございます。 ハンドブックは量が多そうですが、 頑張って調べてみます。