大学の課題に行き詰って…

[ケンイチ]

大学の課題で以下のサイトの人の論文に出てくる数値計算プログラムを書いています。
http://ci.nii.ac.jp/els/110002723163.pd ... 816589&cp=

とりあえず書き終わり、実行結果を見ると、望んだ数値が出ませんでした。ここ何日かデバッグをしても発展がなかったので、今回ここを頼ることにしました。些細なことでも指摘していただけると助かります。なお、書いたプログラムは論文の§3の計算方法を使用し、数値は§5.1の初期設定を使用しています。

コード:

#include<stdio.h>
#define _USE_MATH_DEFINES
#include<math.h>

#define N 100


void Initial_z(int n,double *zx,double *zy,double p,double z_0){//zの生成
	int i;
	for(i=0;i<n;i++){
		zx[i]=p*cos(2*M_PI*(i+1)/n)+z_0;
		zy[i]=p*sin(2*M_PI*(i+1)/n);
	}
}

void Initial_xi(int n,double *xi_x,double *xi_y,double zx[N],double zy[N],double p,double q,double z_0){//ξの生成
	int i;
	for(i=0;i<n;i++){
		if(z_0<zx[i])
			xi_x[i]=zx[i]+fabs((q-p)*cos(2*M_PI*(i+1)/n));
		else
			xi_x[i]=zx[i]-fabs((q-p)*cos(2*M_PI*(i+1)/n));
		xi_y[i]=q*zy[i];
	}
}

void Initial_GN(int n,double GN[N][N],double zx[N],double zy[N],double xi_x[N],double xi_y[N]){//係数行列の生成
	int i,j;
	for(i=0;i<n;i++){
		for(j=0;j<n;j++){
			GN[i][j]=sqrt((zx[i]-xi_x[j])*(zx[i]-xi_x[j])+(zy[i]-xi_y[j])*(zy[i]-xi_y[j]));
			GN[i][j]=log(GN[i][j]);
		}
	}
}

void Initial_logz(int n,double zx[N],double zy[N],double *logz){//連立1次方程式の右辺を生成
	int i;
	for(i=0;i<n;i++)
		logz[i]=log(sqrt(zx[i]*zx[i]+zy[i]*zy[i]));
}

void Gauss_Method(int n,double G[N][N],double f[N],double *Q){//連立一次方程式を解く関数
	int i,j,k;
	double s;	
	/*前進消去*/
	for(i=0;i<n-1;i++){
		for(j=i+1;j<n;j++){
			s=G[j][i]/G[i][i];
			for(k=i+1;k<n;k++)
				G[j][k]=G[j][k]-G[i][k]*s;
			f[j]=f[j]-f[i]*s;
		}
	}
	/*後進消去*/
	for(i=n-1;i>=0;i--){
		s=f[i];
		for(j=i+1;j<n;j++)
			s=s-G[i][j]*f[j];
		f[i]=s/G[i][i];
		Q[i]=f[i];
	}
}
void Calc_G(int n,double Q[N],double zx[N],double zy[N],double xi_x[N],double xi_y[N],double *G){//Gの計算
	int i,j;
	double sum,num;
	for(i=0;i<n;i++){
		sum=0;
		for(j=0;j<n;j++){
			num=sqrt((zx[i]-xi_x[j])*(zx[i]-xi_x[j])+(zy[i]-xi_y[j])*(zy[i]-xi_y[j]));
			sum+=Q[j]*log(num);
		}
		G[i]=-1*sum;
	}
}
void Calc_H(int n,double Q[N],double zx[N],double zy[N],double xi_x[N],double xi_y[N],double *H){//Hの計算
	int i,j;
	double sum,x,y;
	for(i=0;i<n;i++){
		sum=0;
		for(j=0;j<n;j++){
			x=1-(zx[i]*xi_x[j]+zy[i]*xi_y[j])/((xi_x[j]*xi_x[j])+(xi_y[j]*xi_y[j]));
			y=(zx[i]*xi_y[j]-zy[i]*xi_x[j])/((xi_x[j]*xi_x[j])+(xi_y[j]*xi_y[j]));
			sum+=Q[j]*atan2(y,x);
		}
		H[i]=-1*sum;
	}
}

void Calc_UV(int n,double G[N],double H[N],double zx[N],double zy[N],double *U,double *V){//U,Vの計算
	int i;
	for(i=0;i<n;i++){
		U[i]=exp(log(sqrt(zx[i]*zx[i]+zy[i]*zy[i]))+G[i])*cos(atan2(zx[i],zy[i])+H[i]);
		V[i]=exp(log(sqrt(zx[i]*zx[i]+zy[i]*zy[i]))+G[i])*sin(atan2(zx[i],zy[i])+H[i]);
	}
}

int main(){
	int n=16,i;
	double zx[N],zy[N],xi_x[N],xi_y[N],GN[N][N],logz[N],Q[N],u[N],v[N],W[N];//計算上必要な配列
	double G[N],H[N],U[N],V[N];//最終的に求めたいものを格納する配列
	double p=1.0;//拘束点の半径
	double q=1.2;//電荷点の拡大率
	double z_0=-0.25;//正規化条件

	Initial_z(n,zx,zy,p,z_0);//拘束点zの生成
	Initial_xi(n,xi_x,xi_y,zx,zy,p,q,z_0);//電荷点ξの生成
	Initial_GN(n,GN,zx,zy,xi_x,xi_y);//連立1次方程式の係数行列を生成
	Initial_logz(n,zx,zy,logz);//連立1次方程式の右辺を生成
	Gauss_Method(n,GN,logz,Q);//電荷Qを求める。
	Calc_G(n,Q,zx,zy,xi_x,xi_y,G);//Gの計算
	Calc_H(n,Q,zx,zy,xi_x,xi_y,H);//Hの計算
	Calc_UV(n,G,H,zx,zy,U,V);//U,Vの生成
	return 0;
}

[beatle]

とりあえず書き終わり、実行結果を見ると、望んだ数値が出ませんでした。

その実行結果と、望んでいる数値をお示しください。

ここ何日かデバッグをしても発展がなかったので、今回ここを頼ることにしました。

どのようなデバッグ手法を使っていますか？

• printfは使っていますか？
• デバッガは使っていますか？

取り敢えず、様々な場所にprintfを埋め込んで重要な変数の値などを出力してみて、望んだ値になっているかをチェックするというのが簡単なデバッグ方法です。お試しください。
もしそれでも解決が難しいなら、printfを埋め込んだソースコードと、そのソースコードの実行結果と、望んでいる出力の3つを貼り付けてもらえるとさらにアドバイスできると思います。