ソースコードを簡潔に直したいのですが…
有限積分法(参考 http://www.akita-nct.jp/yamamoto/study/thesis/2005/thesis_namekawa.pdf)
を用いて、格子点 ( i, j, k )の磁場成分 Bx, By, Bz をその周辺を囲む電場成分 Ex, Ey, Ez で逐次計算するコードは次のようになります。
for ( i = 0; i < i_max; i++ ) {
for ( j = 0; j < j_max; j++ ) {
for ( k = 0; k < k_max; k++ ) {
Bx [ i ][ j ][ k ] = Ey [ i + 1 ][ j ][ k ] - Ey [ i + 1 ][ j ][ k + 1 ] + Ez [ i + 1 ][ j + 1 ][ k ] - Ez [ i + 1 ][ j ][ k ];
By [ i ][ j ][ k ] = Ez [ i ][ j + 1 ][ k ] - Ez [ i + 1 ][ j + 1 ][ k ] + Ex [ i ][ j + 1 ][ k + 1 ] - Ex [ i ][ j + 1 ][ k ];
Bz [ i ][ j ][ k ] = Ex [ i ][ j ][ k + 1 ] - Ex [ i ][ j + 1 ][ k + 1 ] + Ey [ i + 1 ][ j ][ k + 1 ] - Ey [ i ][ j ][ k + 1 ];
}
}
}
電場成分の係数行列 C ( 0, 1, -1 のいずれかをもつ ) を求める必要が生じた、すなわち
CCC E B
CCC * E = B
CCC E B
のような行列計算に変更するため、次のように書き直しましたが係数行列 C を求める部分を簡潔にできませんでした。
for ( i = 0; i < i_max; i++ ) {
for ( j = 0; j < j_max; j++ ) {
for ( k = 0; k < k_max; k++ ) {
for ( x = 0; x < 3; x++ ) {
row = index ( x, i, j, k, i_max, j_max, k_max );
y = ( x + 1 ) % 3;
z = ( x + 2 ) % 3;
if ( x == 0 ) {
C[ row ][ index ( y, i + 1, j, k, i_max, j_max, k_max )] = 1;
C[ row ][ index ( y, i + 1, j, k + 1, i_max, j_max, k_max )] = -1;
C[ row ][ index ( z, i + 1, j + 1, k, i_max, j_max, k_max )] = 1;
C[ row ][ index ( z, i + 1, j, k, i_max, j_max, k_max )] = -1;
}
if ( x == 1 ) {
C[ row ][ index ( y, i, j + 1, k, i_max, j_max, k_max )] = 1;
C[ row ][ index ( y, i + 1, j + 1, k, i_max, j_max, k_max )] = -1;
C[ row ][ index ( z, i, j + 1, k + 1, i_max, j_max, k_max )] = 1;
C[ row ][ index ( z, i, j + 1, k, i_max, j_max, k_max )] = -1;
}
if ( x == 2 ) {
C[ row ][ index ( y, i, j, k + 1, i_max, j_max, k_max )] = 1;
C[ row ][ index ( y, i, j + 1, k + 1, i_max, j_max, k_max )] = -1;
C[ row ][ index ( z, i + 1, j, k + 1, i_max, j_max, k_max )] = 1;
C[ row ][ index ( z, i, j, k + 1, i_max, j_max, k_max )] = -1;
}
}
}
}
}
int index ( const int d, const int i, const int j, const int k, const int i_max, const int j_max, const int k_max )
{
return ( ( d * i_max + i ) * j_max + j ) * k_max + k;
}
void Solver ( const int i_max, const int j_max, const int k_max)
{
int row, col,
cell_num = i_max * j_max * k_max;
for ( row = 0; row < cell_num; row++ ) {
for ( col = 0; col < cell_num; col++ ) {
B [ row ] += C [ row ][ col ] * E [ col ];
}
}
}
添え字計算が多く、非常に複雑なソースコードですが、規則性があるので、
皆様の力をお借りして簡潔に表現したいです。
よろしくお願いします。
お礼
ご回答ありがとうございます。 >定式化上、整然とするという点がメリットに感じます。 おっしゃる通りです。 >係数行列は、書かれた式をそのまま行列形式で書き下すだけです。 各軸の計算格子点を独立に(空間を立方体でなく直方体に)変更したり、 グリッド形状を変更した際に係数行列が生成できたらと思いまして、 コーディングも苦手ですが頑張ってみます。
補足
「周波数領域で解く場合」とはどのような場合ですか。 放射された電磁波のスペクトルを見る場合などでしょうか。 浅学で申し訳ございませんが、 「周波数領域で解く場合の行列を用いる必然性」を教えて頂けませんか。