Q&A
説明不足でなにを聞きたいのかわかりません。
もちっと詳しく、具体的に説明できませんか
前提
あるk_xに対して固有値を計算するプログラムを書いています。
実現したいこと
固有値自体は計算できているようですが、k_xの変化に対応できていない
「あるk_xに対してこの値」みたいなプログラミングにしたい
発生している問題・エラーメッセージ
同じ値しか出てこない
あるk_xに対して行列a[][]が決まって、そのa[][]の固有値を計算してその値をD[]に入れるようなプログラミングが書きたいのですが、ほしい値が出てきていないのが現状です
該当のソースコード
C
1//ヤコビ法による固有値と固有ベクトル計算// 2 3#include <stdio.h> 4#include <math.h> 5#include <complex.h> 6#include <stdlib.h> 7 8#define N 4//次数設定 9#define L 8 10#define EPS 0.0001 //収束範囲 11 12int n = 100000; 13int M = 2; //電子数＝2M 14int c = 1; 15double complex t_1, t_2, t_3, t_4, t_5, t_6; //方向 16double b_pi; 17double a[L][L]; 18double D[2*N]; 19 20 21 22int main(int argc, char** argv) { 23 24 b_pi = M_PI; 25 double complex a_cube = 0.246 * pow(10, -9); //格子定数[nm] 26 double complex k_x = -b_pi / a_cube; //↓のk_xの値を変えるのも忘れずに！ 27 double complex k_y = b_pi / M * c; 28 29 FILE* fp; 30 31 if ((fp = fopen("zig N=2 q=1.txt", "w")) == NULL) 32 { 33 printf("Cannot open the file\n"); 34 exit(1); 35 } 36 37 int v; 38 for (v = 0;v < n; v++) { 39 k_x = k_x + b_pi * 2 / a_cube / n; 40 t_4 = cos(-1 * a_cube / pow(3, 0.5) * k_y) + I * sin(-1 * a_cube / pow(3, 0.5) * k_y); //cexp(I * (-1* a_cube / pow(3, 0.5) * k_y)); 41 t_5 = cos((a_cube / 2 / pow(3, 0.5) * k_x + a_cube / 2 * k_y)) + I * sin((a_cube / 2 / pow(3, 0.5) * k_x + a_cube / 2 * k_y)); //cexp(I * (a_cube / 2 / pow(3, 0.5) * k_x + a_cube / 2 * k_y)); //expの中の符号を計算したあらかじめ 42 t_6 = cos(a_cube / 2 / pow(3, 0.5) * k_x - a_cube / 2 * k_y) + I * sin(a_cube / 2 / pow(3, 0.5) * k_x - a_cube / 2 * k_y); //cexp(I * (a_cube / 2 / pow(3, 0.5) * k_x - a_cube / 2 * k_y)); // 43 t_1 = cos(a_cube / pow(3, 0.5) * k_y) + I * sin(a_cube / pow(3, 0.5) * k_y); //cexp(I * a_cube / pow(3, 0.5) * k_y); 44 t_2 = cos(-a_cube / 2 / pow(3, 0.5) * k_x - a_cube / 2 * k_y) + I * sin(-a_cube / 2 / pow(3, 0.5) * k_x - a_cube / 2 * k_y); //cexp(I * (-a_cube / 2 / pow(3, 0.5) * k_x - a_cube / 2 * k_y)); 45 t_3 = cos(-a_cube / 2 / pow(3, 0.5) * k_x + a_cube / 2 * k_y) + I * sin(-a_cube / 2 / pow(3, 0.5) * k_x + a_cube / 2 * k_y); //cexp(I * (-a_cube / 2 / pow(3, 0.5) * k_x + a_cube / 2 * k_y)); 46 47 double complex C[N][N] = { 48 {0,0,t_5 + t_6,0}, 49 {0,0,t_4,t_5 + t_6}, 50 {t_2 + t_3,t_1, 0, 0}, 51 {0,t_2 + t_3, 0, 0} }; //係数行列 52 53 double A[N][N], B_1[N][N], B_2[N][N]; 54 int k, h; 55 for (k = 0;k < N;k++) { 56 for (h = 0;h < N;h++) { 57 A[k][h] = creal(C[k][h]); 58 B_1[k][h] = -1 * cimag(C[k][h]); 59 B_2[k][h] = cimag(C[k][h]); 60 } 61 } 62 /*for (k = 0;k < N;k++) { 63 for (h = 0;h < N;h++) { 64 printf("%7.4lf\n", A[k][h]); 65 } 66 }*/ 67 68 69 /*= { 70 {A[0][0],A[0][1],A[0][2] ,B_1[0][0],B_1[0][1],B_1[0][2]}, 71 {A[1][0],A[1][1],A[1][2] ,B_1[1][0],B_1[1][1],B_1[1][2]}, 72 {A[2][0],A[2][1],A[2][2] ,B_1[2][0],B_1[2][1],B_1[2][2]}, 73 {B_2[0][0],B_2[0][1],B_2[0][2] ,A[0][0],A[0][1],A[0][2]}, 74 {B_2[1][0],B_2[1][1],B_2[1][2] ,A[1][0],A[1][1],A[1][2]}, 75 {B_2[2][0],B_2[2][1],B_2[2][2] ,A[2][0],A[2][1],A[2][2]}};*/ 76 for (k = 0;k < L;k++) { 77 for (h = 0;h < L;h++) { 78 if (h < N){ 79 if (k < N) 80 a[h][k] = A[h % N][k % N]; 81 else 82 a[h][k] = B_1[h % N][k % N]; 83 } 84 else { 85 if (k < N) 86 a[h][k] = B_2[h % N][k % N]; 87 else 88 a[h][k] = A[h % N][k % N]; 89 } 90 } 91 } 92 93 /*/for (k = 0;k < L;k++) { 94 for (h = 0;h < L;h++) 95 printf("%7.4lf\n", a[k][h]); 96 }*/ 97 98 //double u[L][L]; //単位行列 99 double alpha, beta, gamma; 100 double s, c, w; 101 double wa, wb, wc; 102 double max; 103 int i, j, p, q, x, y; 104 105 /*for (y = 0;y < N;y++) 106 for (x = 0;x < N;x++) 107 u[y][x] = 0.0; 108 109 for (i = 0;i < N;i++) 110 u[i][i] = 1.0;*/ 111 112 while (1) { 113 //最大要素の行と列を検索 114 max = 0.0; 115 for (i = 0;i < L - 1;i++) { 116 for (j = i + 1;j < L;j++) 117 if (fabs(a[i][j]) > max) { 118 p = i; 119 q = j; 120 max = fabs(a[i][j]); 121 } 122 } 123 //収束したら解答打ち出し 124 if (max < EPS) break; 125 126 //sin, cos 計算 127 wa = a[p][p]; 128 wb = a[p][q]; 129 wc = a[q][q]; 130 alpha = -wb; 131 beta = 0.5 * (wa - wc); 132 gamma = fabs(beta) / sqrt(alpha * alpha + beta * beta); 133 s = sqrt(0.5 * (1.0 - gamma)); 134 if (alpha * beta < 0) s = -s; 135 c = sqrt(1.0 - s * s); 136 //直行変換 137 for (j = 0;j < L;j++) { 138 w = a[p][j] * c - a[q][j] * s; 139 a[q][j] = a[p][j] * s + a[q][j] * c; 140 a[p][j] = w; 141 } 142 143 for (j = 0;j < L;j++) { 144 a[j][p] = a[p][j]; 145 a[j][q] = a[q][j]; 146 147 148 } 149 150 w = 2.0 * wb * s * c; 151 a[p][p] = wa * c * c + wc * s * s - w; 152 a[q][q] = wa * s * s + wc * c * c + w; 153 a[p][q] = 0; 154 a[q][p] = 0; 155 156 //漸化式計算 157 /*for (i = 0;i < N;i++) { 158 w = u[i][p] * c - u[i][q] * s; 159 u[i][q] = u[i][p] * s + u[i][q] * c; 160 u[i][p] = w; 161 }*/ 162 163 D[0] = a[0][0]; 164 D[1] = a[1][1]; 165 D[2] = a[2][2]; 166 D[3] = a[3][3]; 167 D[4] = a[4][4]; 168 D[5] = a[5][5]; 169 D[6] = a[6][6]; 170 D[7] = a[7][7]; 171 172 173 /*for (k = 0;k < L;k++) { 174 for (h = 0;h < L; h++) { 175 printf("%7.4lf\n", a[k]); 176 } 177 }*/ 178 } 179 180 fprintf(fp, "%2f %lf %lf %lf %lf %lf %lf %lf %lf\n", k_x, D[0], D[1], D[2], D[3], D[4], D[5], D[6], D[7]); 181 } 182 183 /*int i; 184 printf("固有値\n"); 185 for (i = 0;i < L;i++) { 186 printf("%7.4lf\n", a[i][i]); //対角成分 こいつが固有値 187 }*/ 188 189 fclose(fp); 190 return 0; 191}
試したこと
ここに問題に対して試したことを記載してください。
補足情報（FW/ツールのバージョンなど）
ここにより詳細な情報を記載してください。
すいません。
あるk_xに対して行列a[][]が決まって、そのa[][]の固有値を計算してその値をD[]に入れるようなプログラミングが書きたいのですが、ほしい値が出てきていないのが現状です
コードのコメントでは最後に固有値を表示しているらしい個所はありますが、そこで D[] に入れれば良いだけでは。
なぜ a の計算途中？のループの中で D[] に入れているのでしょうか。(表示はループの外ですし。)
(なんとか法とか全く知らないので何処で何をやるのが正解なのか知りませんが。)
すいません、今のコードを書きます。
自己解決
