固有値の数値計算 複素数を含め

C

//ヤコビ法による固有値と固有ベクトル計算// #include <stdio.h>#include <math.h>#include <stdlib.h>#include <complex.h> #define N 4 //次数設定#define EPS 0.0001 //収束範囲 int k,r,t;int n=100000;int M=2; //電子数＝2Mint c=1; //波数 c=1,2,...2Mdouble a_cube=0.246*pow(10,-9); //格子定数[nm]double b_pi=M_PI; //円周率 int main( int argc, char **argv) { double k_x=-b_pi/a_cube; //範囲を拡張するのには/2を消せばいいじゃない。↓のk_xの値を変えるのも忘れずに！ double k_y=b_pi/M*c; FILE *fp; if ((fp=fopen("zig N=2.txt","w")) == NULL){ printf("Cannot open the file\n"); exit(1); } for(k=0;k<100000;k++){ k_x=k_x+b_pi*2/a_cube/n; double a[N][N]={ { 0,0,exp(-I*(-a_cube/2/pow(3,0.5)*k_x-a_cube/2*k_y))+exp(-I*(-a_cube/2/pow(3,0.5)*k_x+a_cube/2*k_y)),0}, {0,0,exp(-I*a_cube/pow(3,0.5)*k_y),exp(-I*(-a_cube/2/pow(3,0.5)*k_x-a_cube/2*k_y))+exp(-I*(-a_cube/2/pow(3,0.5)*k_x+a_cube/2*k_y))}, {exp(I*(-a_cube/2/pow(3,0.5)*k_x-a_cube/2*k_y))+exp(I*(-a_cube/2/pow(3,0.5)*k_x+a_cube/2*k_y)),exp(-I*a_cube/pow(3,0.5)*k_y),0,0}, {0,exp(I*(-a_cube/2/pow(3,0.5)*k_x-a_cube/2*k_y))+exp(I*(-a_cube/2/pow(3,0.5)*k_x+a_cube/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[0][0],a[1][1],a[2][2],a[3][3]); fprintf(fp, "%2f %lf %lf %lf %lf\n",k_x,a[0][0],a[1][1],a[2][2],a[3][3]); } } fclose(fp); return 0;}