ニュートン法のプログラムについてです。

[ごんた]

こんばんは。
今、ニュートン法のプログラムをやっているのですが、僕の力では解決できないエラーがでて途方に暮れている状態です。
恐らく、ガウスの消去法を行っているところでエラーがでていると思うのですが、どこをどう直せばよいかわかりません。
どなたかわかる方、教えてくださると助かります。
参考にしたサイトはhttp://pc-physics.com/gauss1.htmlです。
以下、プログラムです。それでは失礼します。

#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#define N 2 /*行列のサイズ*/

/*超越関数 f(x,y)=x+y^2*/
double f(double x, double y)
{
return x+y*y;
}

/*f一回x偏微分*/
double fdx(double x,double y)
{
return 1;
}

/*f一回y偏微分*/
double fdy(double x, double y)
{
return 2*y;
}

/*超越関数 g(x,y)=x^2+(y-1)^2-1*/
double g(double x, double y)
{
double z;

z=x*x+(y-1)*(y-1)-1;
return z;
}

/*g一回x偏微分*/
double gdx(double x, double y)
{
return 2*x;
}

/*g一回y偏微分*/
double gdy(double x, double y)
{
return 2*(y-1);
}


/*メインプログラム*/
int main()
{
int i,j,k,l,pivot;
double a[N][N+1];/*行列a[行][列]*/
double h[N];
double x=2.57;
double y=3.43;
double buf,p,m,b[1][N+1];;

h[0]=0;
h[1]=0;

while((fabs(f(x,y))>0.0000001)&&(fabs(g(x,y))>0.0000001))/*f,gがともに一
定数より大きかったらwhileループ*/
{

a[0][0]=fdx(x,y);
a[0][1]=fdy(x,y);
a[1][0]=gdx(x,y);
a[1][1]=gdy(x,y);
a[0][2]=f(x,y);
a[1][2]=g(x,y);

for(i=0;i<N;i++){

m=0;

pivot=i;

for(l=i;l<N;l++){

if(fabs(a[l])>m){

m=fabs(a[l]);

pivot=l;

}
}

if(pivot!=i){
for(j=0;j<N+1;j++){

b[0][j]=a[j];
a[j]=a[pivot][j];
a[pivot][j]=b[0][j];

}
}

}



/*左下三角を0にしたい。*/
for(i=0;i<N;i++){

buf=1/a; /*対角要素を保存*/

a=1; /*対角要素は1になるので直で代入。*/

for(j=i+1;j<N;j++){

a[j]=a[j]*buf; /*同じ行の成分にかけまくり*/

}

for(k=i+1;k<N+1;k++){

p=a[k][i];

for(j=i+1;j<N+1;j++){

a[k][j]=a[k][j]-p*a[i][j];
}

a[k][i]=0; /*0になることがわかってるので直で代入*/

}
}/*行列の諸々おわり*/

/*解の計算*/

for(i=N-1;i>=0;i--){

h[i]=a[i][N]; /*階段の最初のところ。*/

for(j=N-1;j>i;j--){

h[i]=h[i]-a[i][j]*h[j];

}
}

x=x+h[0];
y=y+h[1];

}

/*結果の表示*/
printf("x=%10.8lf\n",x);
printf("y=%10.8lf\n",y);
printf("f(x,y)=%10.8lf\n",f(x,y));
printf("g(x,y)=%10.8lf\n",g(x,y));
}

[box]

今、ニュートン法のプログラムをやっている

解きたい問題の中身や、出ているエラーの具体的な内容に関する説明は、ナシですか？
その説明もナシにいきなりコードを見ろ、といわれましてもね…。

何を書けば回答が得られそうか、事前に想定しましょう。

[bitter_fox]

ごんたさん、フォーラムルールはお読みになられましたか？
ソースコードを乗せる際はcodeタグで囲むようにとあります。これに従うとインデントが再現され非常に読みやすくなります。ですのでこれに従っていただくようにお願いします。

また具体的にどういったエラーが出ているのかを教えていただくと回答者の方々が答えやすくなりますのでもう少し詳しく教えていただけますか？

コード:

#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#define N 2 /*行列のサイズ*/

/*超越関数 f(x,y)=x+y^2*/
double f(double x, double y)
{
	return x+y*y;
}

/*f一回x偏微分*/
double fdx(double x,double y)
{
	return 1;
}

/*f一回y偏微分*/
double fdy(double x, double y)
{
	return 2*y;
}

/*超越関数 g(x,y)=x^2+(y-1)^2-1*/
double g(double x, double y)
{
	double z;

	z=x*x+(y-1)*(y-1)-1;
	return z;
}

/*g一回x偏微分*/
double gdx(double x, double y)
{
	return 2*x;
}

/*g一回y偏微分*/
double gdy(double x, double y)
{
	return 2*(y-1);
}


/*メインプログラム*/
int main()
{
	int i,j,k,l,pivot;
	double a[N][N+1];/*行列a[行][列]*/
	double h[N];
	double x=2.57;
	double y=3.43;
	double buf,p,m,b[1][N+1];;

	h[0]=0;
	h[1]=0;

	while((fabs(f(x,y))>0.0000001)&&(fabs(g(x,y))>0.0000001))/*f,gがともに一
	定数より大きかったらwhileループ*/
	{

		a[0][0]=fdx(x,y);
		a[0][1]=fdy(x,y);
		a[1][0]=gdx(x,y);
		a[1][1]=gdy(x,y);
		a[0][2]=f(x,y);
		a[1][2]=g(x,y);

		for(i=0;i<N;i++){

			m=0;

			pivot=i;

			for(l=i;l<N;l++){

				if(fabs(a[l][i])>m){

					m=fabs(a[l][i]);

					pivot=l;

				}
			}

			if(pivot!=i){
				for(j=0;j<N+1;j++){

					b[0][j]=a[i][j];
					a[i][j]=a[pivot][j];
					a[pivot][j]=b[0][j];

				}
			}

		}



	/*左下三角を0にしたい。*/
		for(i=0;i<N;i++){

			buf=1/a[i][i]; /*対角要素を保存*/

			a[i][i]=1; /*対角要素は1になるので直で代入。*/

			for(j=i+1;j<N;j++){

				a[i][j]=a[i][j]*buf; /*同じ行の成分にかけまくり*/

			}

			for(k=i+1;k<N+1;k++){

				p=a[k][i];

				for(j=i+1;j<N+1;j++){

				a[k][j]=a[k][j]-p*a[i][j];
				}

				a[k][i]=0; /*0になることがわかってるので直で代入*/

			}
		}/*行列の諸々おわり*/

		/*解の計算*/

		for(i=N-1;i>=0;i--){

			h[i]=a[i][N]; /*階段の最初のところ。*/

			for(j=N-1;j>i;j--){

				h[i]=h[i]-a[i][j]*h[j];

			}
		}

		x=x+h[0];
		y=y+h[1];

	}

	/*結果の表示*/
	printf("x=%10.8lf\n",x);
	printf("y=%10.8lf\n",y);
	printf("f(x,y)=%10.8lf\n",f(x,y));
	printf("g(x,y)=%10.8lf\n",g(x,y));
}

[box]

まあ、ザッと見ただけですけどね。

コード:

    /*左下三角を0にしたい。*/
        for(i=0;i<N;i++){
            buf=1/a[i][i]; /*対角要素を保存*/
            a[i][i]=1; /*対角要素は1になるので直で代入。*/
            for(j=i+1;j<N;j++){
                a[i][j]=a[i][j]*buf; /*同じ行の成分にかけまくり*/
            }
            for(k=i+1;k<N+1;k++){
                p=a[k][i];
                for(j=i+1;j<N+1;j++){
                a[k][j]=a[k][j]-p*a[i][j];
                }
                a[k][i]=0; /*0になることがわかってるので直で代入*/

外側の、iに関するループで、iはN-1の値を持つことがあるので、
内側の、kに関するループで、kの初期値がi+1ってことは、
kがNという値を持つことがある、ってことですよね。
配列aの定義範囲を超えてませんか？

[ごんた]

申しわけありません。おっしゃるとおりです。
フォーラムルールも読んでいませんでした。ご迷惑をおかけして本当に申し訳ありません。

解きたい内容は連立方程式
x+y~2=0
x~2+(y-1)~2-1=0
をニュートン法のアルゴリズムを用いて解を求めよというものです。
エラーは、コンパイルして実行した後、

Debug Error!

Program: ...（ファイル名）
Module: ...(ファイル名)
File:

Run-Time Check Failure #2 - Stack around the variable 'a' was corrupted.

(Press Retry to debug the application)

という警告(?)ウインドウが出てきてしまいます。
また実行結果も最後のところが
x=-1.#IND0000
y=-1.#IND0000
f(x,y)=-1.#IND0000
g(x,y)=-1.#IND0000
となり、うまくいきません。
0や小さい数などで除算を行っているときなどにこのような結果になるらしいのですが、プログラムのどこの段階でそのようなことをやっているかもわからない状態です。
OSはWindows Vistaで、コンパイラはVC2010Expressです。

[ごんた]

box様
レスありがとうございます。
ご指摘いただいたところを修正したところ、警告ウインドウはでなくなりました。
ただ、結果がやはりおかしくなります。

[たいちう]

どうおかしくなるのか書かないと。
どういう結果を期待しているのかと、どう修正したのかも。

[box]

ただ、結果がやはりおかしくなります。

こういう風に書くから、「それで、どんな風に？」と言われるのです。
具体的に、どうおかしいかを一発で書いてくだされば、こういうよけいなやりとりをしなくてすむんですけどね。

[ごんた]

解決しました。

コード:

for(j=i+1;j<N;j++){
 
                a[i][j]=a[i][j]*buf; /*同じ行の成分にかけまくり*/
 
            }

のjの条件で、
j<N⇒j<N+1
としたら上手くいきました。
お手数をおかけして申し訳ありませんでした。
box様、bitter_fox様、たいちう様、お忙しいところご助言をくださいましてありがとうございました。

たいちう様、box様
レスありがとうございます。
結果がおかしくなるというのは、その前に僕が書き込みした
x=...
y=...
と同じだったので、省略したつもりで書いてしまいました。わかりづらくて申しわけありません。