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;
}

回答No.2

noname#88356 · ベストアンサー

double *simple_gauss( double a[3][3], double b[3] );
などはこれでも正解なのですが、具体的な数値名が関数の定義に含んでいるので
本当はポインタであることを明確にするために
double *simple_gauss( double **a, double *b );
などとするべきかなと思います。（その方が一々数字を書き換えなくていいですし）

後、細かいことではありますが、
double a[3][3]; /* マトリックスの宣言*/
double x[3]; /* 解ベクトルの宣言*/
double y[3]; /* 右辺ベクトルの宣言*/
double xx, yy; /* キーボード入力のための一時変数*/

の部分は個々の変数名が分かりにくいのも問題かなと思います。
私も昔はこういう書き方をしていましたが、エラーが出る可能性の増大が馬鹿に出来ないことに
ある時気がつき、forループなどの一時的なカウンター変数を除いては一時的なものであっても
分かりやすい命名を心がけています。

さらに、これらの変数は見たところ共通している要素であるため、構造体、または
Ｃ＋＋であればクラスなどでパックにした方が後々の修正や読む人が楽であると思います。

最後に、コメントに
/* コメント */
の方法を取っているようですが、これは単一の行で用いるにはあまり推奨されない方法です。
というのも/*　*/はその間に/* */をはさんでコメントアウトなどが出来ないので、
もしある関数の中身でコメントに/*　*/が多用されているとテストするときに
ここからここまでコメントアウトしたい・・などというときに/*　*/を使って一気に
コメントアウトしようとしても間に/*　*/が混じっているとコンパイルエラーになってしまい、
大規模なソフトを作る際にはかなり面倒が起こりやすいです。
初期のＣ言語では/*　*/が唯一のコメント方法であり、//コメント の形式はサポートされていませんでしたが、
現在は標準Ｃで//コメントの形式がサポートされているので、例えば

double a[3][3]; /* マトリックスの宣言*/
double x[3]; /* 解ベクトルの宣言*/
double y[3]; /* 右辺ベクトルの宣言*/
double xx, yy; /* キーボード入力のための一時変数*/

は

double a[3][3]; //マトリックスの宣言
double x[3];　　//解ベクトルの宣言
double y[3];　　//右辺ベクトルの宣言
double xx, yy;　//キーボード入力のための一時変数

のようにした方が良いかと思います。

後、要素の数が３の配列が多いですが、これらは必ず３と固定されているのでしょうか？
もし数学的な理由で３と固定されているなら

#define ARY_NUM 3
などとして配列の宣言時は
double a[ ARY_NUM ][ ARY_NUM ]; 　//マトリックスの宣言
double x[ ARY_NUM ];　　　　　　　　　　//解ベクトルの宣言
double y[ ARY_NUM ];　　　　　　　　　 //右辺ベクトルの宣言
double xx, yy;　　　　　　　　　　　　　　　//キーボード入力のための一時変数
などとした方が色々な面で後の変更や修正が楽です
（Ｃ＋＋ならconst intでやる方法もあります）

ところで自宅のＰＣにはVisual Studioを入れていないということでしたが、
http://www.microsoft.com/japan/msdn/vstudio/express/
無料版があることは既にご存知でしょうか？
おせっかいかもしれませんが、もしご存知でなければ本当に便利ですので導入をオススメします。
（私もCやC++で３Ｄのプログラム組む時などはこれが手放せないです）

お礼

Visual Studioの情報ありがとうございます

aが3×3の行列なので、縦横の順にa[3][3]とおきました
対応してx、yは3×1行列なのでx[3]、y[3]とおきました
a×x＝yの関係です

求めるa,b,cと行列aが同じに見えてややこしいのでp,q,rに変えますね、すいません
要は以下の連立方程式を行列化しただけなのですが

p×Σxi^4＋q×Σxi^3＋r×Σxi^2＝Σxi^2*yi
p×Σxi^3＋q×Σxi^2＋r×Σxi＝Σxi*yi
p×Σxi^2＋q×Σxi＋r×5＝Σyi
見ての通りこれを行列化しa×x＝yにするとaがややこしくなるのでa[1][1]＝～,a[1][2]＝～....と成分を書き出しました

main(void)内で5点の座標(x,y)をキーボード入力して、その後の計算して表示double *simple_gaussで計算した結果を表示してくれればいいのですが

それで
http://www.ma.is.saga-u.ac.jp/minamoto/book/book8/book8.html
の本を授業で使ってまして、これの第３章3.1のプログラムをもとに作成したのですが、これはinput.dat等のファイルから展開して計算をしているのですがこれをscanfで直接5点の座標(x,y)をキーボード入力して取り込みそれをもとに計算するように変えようとして失敗したのが始めに投稿したプログラムです
ガウス消去法の部分はそのまま使えるのでmain(void)部分内を何とか変えることができたらと思うんですが

ちなみに/* コメント */もこの本を参考にしているので詳しいことは知りませんでした、全くの素人なんでfatalitaさんのような専門の人に質問を回答いただいて本当に感謝しています

回答No.1

noname#88356

私には数学の知識が中卒程度しかないのでプログラム部分だけアドバイスを

・main(void)はint main( void )と書くべき
・simple_gauss関数でfor( i = k+1; i <= ; i++)となっている（タイプミスでは？
・ y=*simple_gauss(a,y);は定義済みのyに配列確保しようとしています。
　mallocするかあるいはdouble * y2 = simple_gauss( a, y );のように新しい配列を確保しなおす必要があります。

これを修正したらコンパイルは通りましたが、これはあまり良いコードとは思えません。

お礼

スペースを空けたのに何故か全部左詰めになって見にくくてすみません。
あまりよいコードではないとは？
点(x,y)を５点キーボード入力してa,b,cがでればいいのですが....

すんません、今返信してるパソコンにはWindowsVisualStudio等のコンパイラがないので大学のパソコンでしか動作確認ができない状態なんですがOTLすいません。

補足

すみませ～ん！ガウス消去法の部分がちょっと間違ってました

{
　　　int i, j, k;
　　　double alpha, tmp;

　　　/* 前進消去 */
　　　for( k = 1; k <= 2; k++)
　　　{
　　　　　for( i = k+1; i <= 3 ; i++)
　　　　　{
　　　　　　alpha = - a[i][k]/a[k][k];
　　　　　　for( j = k+1; j <= 3; j++)
　　　　　　{
　　　　　　　a[i][j] = a[i][j] + alpha * a[k][j];
　　　　　　}
　　　　　　y[i] = y[i] + alpha * y[k];
　　　　　}
　　　　}

　　　/* 後退代入 */
　　　y[3] = y[3]/a[3][3];
　　　for( k = 2; k >= 1; k--)
　　　{
　　　　　tmp = y[k];
　　　　　for( j = k+1; j <= 3; j++)
　　　　　{
　　　　　　tmp = tmp - a[k][j] * y[j];
　　　　　}
　　　　y[k] = tmp/a[k][k];
　　　}
　　　return y;
}

です

仰るとおりタイプミスも修正しましたが....