複素数を含めて固有値計算がしたい 数値計算

Question

前提

オイラー方程式を含んだ行列の固有値の計算

実現したいこと

エネルギーバンド図を書くために固有値を計算したい

発生している問題・エラーメッセージ

多すぎてちょっとわからないです。

該当のソースコード

//ヤコビ法による固有値と固有ベクトル計算//

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <complex.h>

#define N  4		//次数設定
#define EPS 0.0001  	//収束範囲

int main( int argc, char **argv) {
   	int k;
    int n=100000;
	int M=2;  //電子数＝2M
	int c=1;   //波数 c=1,2,...2M
	double a=0.246*pow(10,-9); //格子定数[nm]
	double b=3.1415926; //円周率
	double k_x=-b/2/a;	//範囲を拡張するのには/2を消せばいいじゃない。↓のk_xの値を変えるのも忘れずに！
	double k_y=b/M*c; 
	double t=2.75;  //飛び移り積分
	double E_1,E_2;  //エネルギー ルートのプラスマイナスで区別
	double e_1,e_2;
	FILE *fp;


	if ((fp=fopen("zig N=2.txt","w")) == NULL){
		printf("Cannot open the file\n");
		exit(1);
	}

   for(k=0;k<n;k++){
	double  a[N][N]={
		    	{ 0,0,exp(-I*(-a/2/pow(3,0.5)*k_x-a/2*k_y))+exp(-I*(-a/2/pow(3,0.5)*k_x+a/2*k_y)),0},
		    	{0,0,exp(-I*a/pow(3,0.5)*k_y),exp(-I*(-a/2/pow(3,0.5)*k_x-a/2*k_y))+exp(-I*(-a/2/pow(3,0.5)*k_x+a/2*k_y))},
		    	{exp(I*(-a/2/pow(3,0.5)*k_x-a/2*k_y))+exp(I*(-a/2/pow(3,0.5)*k_x+a/2*k_y)),exp(-I*a/pow(3,0.5)*k_y),0,0},
			{0,exp(I*(-a/2/pow(3,0.5)*k_x-a/2*k_y))+exp(I*(-a/2/pow(3,0.5)*k_x+a/2*k_y)),0,0}};		//係数行列
	double u[N][N];					//単位行列
	double alpha, beta, gamma;
	double s,c,w;
	double wa,wb,wc;
	double max;
	int i, j, p, q, x, y;

	    for(y=0;y<N;y++)
	    	for(x=0;x<N;x++)
			    u[y][x]=0.0;
	
    	for(i=0;i<N;i++)
	    	u[i][i]=1.0;

    	while(1){
								//最大要素の行と列を検索                                                       
	    	max=0.0;
	    	for(i=0;i<N-1;i++)
	    		for(j=i+1;j<N;j++)
	    			if(fabs(a[i][j])>max){
	    				p=i;
	    				q=j;
	    				max=fabs(a[i][j]);
	    			}
								//収束したら解答打ち出し
	    	if(max<EPS) break;
			
								//sin, cos　計算
	    	wa=a[p][p];
	    	wb=a[p][q];
		    wc=a[q][q];
    		alpha=-wb;
    		beta=0.5*(wa-wc);
	    	gamma=fabs(beta)/sqrt(alpha * alpha + beta * beta);
    		s=sqrt(0.5*(1.0-gamma));
    		if(alpha*beta<0) s=-s;
    		c=sqrt(1.0-s*s);
								//直行変換
    		for(j=0;j<N;j++){
    			w=a[p][j]*c-a[q][j]*s;
    			a[q][j]=a[p][j]*s+a[q][j]*c;
    			a[p][j]=w;
    		}

    		for(j=0;j<N;j++){
    			a[j][p]=a[p][j];
    			a[j][q]=a[q][j];


    		}

    		w=2.0*wb*s*c;
    		a[p][p]=wa*c*c+wc*s*s-w;
    		a[q][q]=wa*s*s+wc*c*c+w;
    		a[p][q]=0;
    		a[q][p]=0;
			
								//漸化式計算
    		for(i=0;i<N;i++){
    			w=u[i][p]*c-u[i][q]*s;
    			u[i][q]=u[i][p]*s+u[i][q]*c;
    			u[i][p]=w;
    		}
	printf("%2f  %lf %lf %lf %lf\n",k_x,a[1],a[2],a[3],a[4]);
	fprintf(fp, "%2f  %lf %lf %lf %lf\n",k_x,a[1],a[2],a[3],a[4]);
	        }
 }
	
	
	printf("固有値\n");
	for(i=0;i<N;i++)
		printf("%7.4lf", a[i][i]);			//対角成分　こいつが固有値
	printf("\n\n");
	printf("固有ベクトル\n");
	for(i=0;i<N;i++){
		for(j=0;j<N;i++)
			printf( "%7.4lf", u[i][j]);		//固有ベクトル
		printf("\n");
    }
	fclose(fp);
	return 0;
}

試したこと

補足情報（FW/ツールのバージョンなど）

グラフェンエネルギーバンド

Comments

[jimbe]

大量のメッセージが出るとビックリするなり悲しくなるなりは理解しますが、やらなければならないことなら、行動するしかありません。
エラーメッセージの最初の 10 件程度でも提示してください。

また、全ての処理を main 内に書くのではなく、処理毎に関数にする、関連データを構造体にする等で処理・データ各々の分割統治を考えた方が良いと思います。

[kisara11235]

コメントありがとうございます。
ちょっと心が折れていました。申し訳ございません。
ある程度解決したので現在、調整しております。

Answer 1

• double aとdouble a[N][N] で名前が被っています。
• 最後の方のfor文で未定義のiが使われています。

---

多すぎてちょっとわからないです。

それでも読みましょう。今回の場合も量は多いですけど内容はほぼ同じことを続けて言われています。

Comments

[kisara11235]

本当に申し訳ないです。ちょっと心がおれていました。

[ozwk]

ちょっと書いてコンパイルが通るかを確かめて、またちょっと書いてコンパイルが通るか確かめて...を繰り返しましょう
私が同じプログラムを書こうと思ったら完成までに100回ぐらいコンパイルすると思います

そうすればいきなりエラーメッセージが1000行出てこようが最後にコンパイルエラーが出てから書き足したところにエラーの原因があるはずですから原因究明が容易になります

いっぺんに書いてたくさんエラーを出してしまった場合でも、
冒頭だけ残してコンパイル通す、通ったらもう少し書き戻してコンパイルを通す...を繰り返して同じようなことができます。