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

[AK]

大学の課題で以下の論文に出てくるアルゴリズムを用いてプログラムを作っています。しかし、実行結果が出てほしい理想値とは全く違う値が出てしまい、その後も試行錯誤してみましたが、行き詰ったので質問します。

http://ci.nii.ac.jp/els/110002724289.pd ... 889389&cp=

以下の数値の左側は私のプログラムにより出た実行結果で、右側は理想値です。
E_M=1.019629e+000 E_M=5.5E-3
E_A=2.471952e-001 E_A=3.0E-3

なお、数値のデバックの箇所はプログラム上に表記しているので、そちらを参照してください。また、このプログラムは複素平面上での計算を行っています。

コード:

#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]));
			printf("%f ",GN[i][j]);//係数行列の表示
			GN[i][j]=log(GN[i][j]);
		}
		printf("\n");
	}
}

void Initial_logz(int n,double zx[N],double zy[N],double *Im_num){//複素数の絶対値を求める
	int i;
	for(i=0;i<n;i++)
		Im_num[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]);
	}
}

void Calc_center_point(int n,double A[N],double *B){//中間点を求める
	int i;
	for(i=0;i<n-1;i++)
		B[i]=(A[i]+A[i+1])/2;
}

double Calc_E_M(int n,double W[N]){//誤差E_Mを求める
	int i;
	double max=0;
	for(i=0;i<n-1;i++){
		if(max < fabs(W[i]-1))
			max=fabs(W[i]-1);
	}
	return max;
}

double Calc_E_A(int n,double H[N]){
	int i;
	double max=0;
	for(i=0;i<n;i++){
		if(max < fabs(H[i]))
			max=fabs(H[i]);
	}
	return max;
}

int main(){
	int i;
	/*-----初期設定-----------------*/
	int n=16;//拘束点・電荷点の数
	double p=1.0;//拘束点の半径
	double q=1.2;//電荷点の拡大率
	double z_0=-0.25;//正規化条件(z_0の数値分、zとξの中心がx軸で移動)
	/*------------------------------*/
	double E_M,E_A;//誤差を表す
	double zx[N],zy[N];//複素平面上での拘束点zの座標
	double xi_x[N],xi_y[N];//電荷点ξの座標
	double Q[N];//電荷Q
	double logz[N];//電荷Qを連立1次方程式で求める際に使う右辺
	double GN[N][N];//電荷Qを求める際に使う係数行列
	double G[N];//調和関数G
	double H[N];//Gの共役調和関数
	double U[N],V[N];//誤差評価で使用
	double u[N],v[N];//U_iとU_{i+1}の中間点,V_iとV_{i+1}の中間点
	double W[N];//W=U+iV;
	Initial_z(n,zx,zy,p,z_0);//拘束点zの生成
	for(i=0;i<n;i++)
		printf("%12f\t%12f\n",zx[i],zy[i]);//拘束点zを表示
	getchar();
	printf("\n");
	Initial_xi(n,xi_x,xi_y,zx,zy,p,q,z_0);//電荷点ξの生成
	for(i=0;i<n;i++)
		printf("%12f\t%12f\n",xi_x[i],xi_y[i]);//電荷点ξを表示
	getchar();
	Initial_GN(n,GN,zx,zy,xi_x,xi_y);//連立1次方程式の係数行列を生成
	Initial_logz(n,zx,zy,logz);//連立1次方程式の右辺を生成
	getchar();
	for(i=0;i<n;i++)
		printf("a%12f\n",logz[i]);//連立一次方程式を表示
	getchar();
	/*GN[0][0]=1,GN[0][1]=-1,GN[0][2]=1;//連立一次方程式を求めるのに使うガウスの消去法の動作チェック用の係数行列
	GN[1][0]=2,GN[1][1]=1,GN[1][2]=-3;   
	GN[2][0]=3,GN[2][1]=2,GN[2][2]=-1;
	logz[0]=-5,logz[1]=19,logz[2]=16;*/
	Gauss_Method(n,GN,logz,Q);//電荷Qを求める。
	for(i=0;i<n;i++)//電荷Qを表示
		printf("%12f\n",Q[i]);
	printf("\n");
	Calc_G(n,Q,zx,zy,xi_x,xi_y,G);//Gの計算
	Calc_H(n,Q,zx,zy,xi_x,xi_y,H);//Hの計算
	for(i=0;i<n;i++)//調和関数G・共役調和関数Hの表示
		printf("%12f\t%12f\n",G[i],H[i]);
	getchar();
	Calc_UV(n,G,H,zx,zy,U,V);//U,Vの生成
	Calc_center_point(n,U,u);//中間点のUを求める
	Calc_center_point(n,V,v);//中間点のVを求める
	Initial_logz(n,u,v,W);//|W|を求める
	E_M=Calc_E_M(n,W);//E_Mを求める。
	E_A=Calc_E_A(n,H);
	printf("E_M=%e\n",E_M);//E_Mの表示
	printf("E_A=%e\n",E_A);//E_Mの表示
	getchar();
	return 0;
}

[beatle]

E_MとE_A以外の各変数の値は理想の値になっているのですか？

[AK]

参照した論文にはE_M,E_A位しか理論値は載ってなかったです。
ただし、拘束点z・電荷点ξの配置はgnuplotで配置を確認。
係数行列GN,連立一次方程式の右辺logzの計算式は対数はサンプルの数値を与えて確認して、手計算と一致したことを確認。
電荷Qを求めるガウスの消去法はサンプルの数値で確認済して、手計算と一致したことを確認。

GN,logz,Qを求める計算式には対数が含まれているので、拘束点・電荷点のデータをそのまま使うと、手計算が追い付かないので、サンプル数値を用いました。

[たいちう]

プログラミングとしては多分難しくないんだけど、
AKさんのお役に立つためには、件の論文を理解して
コーディングを検証する作業が必要ですね。
つまりAKさんが本来やるレベルの質と量の作業。

出来る人もこの掲示板にはざらにいると思うけど、
それなりに時間がかかるかも。
ちなみに締め切りとかはありますか？

[AK]

締め切りは来週の月曜になります。
とりあえず、自分も、もう一度論文を読み返してみます。

[みょん]

もしかしたらもう解決されたかもしれませんが、念の為。
E_M = max ||W(z_{j+1/2})|-1|, z_{j+1/2} \in C
と論文にはありますが、Calc_E_M函数の中身(111行目です)では

コード:

            max=fabs(W[i]-1);

となっており、これはmax|W(z)-1|を表しています
なので、ここを

コード:

            max=fabs(fabs(W[i])-1);

と変えてみたら、
E_M=9.803710e-01
E_A=2.471952e-01
とそれっぽい値になったので、ここのミスではないかと思われますがいかがでしょうか？