プログラマ専用SNS ミクプラ

コード:

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

コード:

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