Jacobi法 をテンプレートにして作成

開始行:
RIGHT:寄稿：東條遼平 * 固有値と固有ベクトルを求める [#d575c2f8] Jacobi法(ヤコビ法)を用いて行列の固有値と固有ベクトルを求... * 固有値、固有ベクトルとは [#rdf1bce0] 任意のn次正方行列Aに対し、 #ref(img53.png,nolink) が成り立つときλを固有値、xを固有ベクトルという。 つまり、あるベクトルxに行列Aを掛けるとベクトルxが定数倍に... * Jacobi法について [#tb97b0dd] #ref(img54.png,nolink) となるとき、行列Aと行列Bの固有値は一致します。また、行列A... #ref(img55.png,nolink) と、対角化することができ、行列Λの対角成分は行列Λの固有値... このことから行列Aの固有値と固有ベクトルを求めるには、行列... Jacobi法ではこの直交行列Pを反復作業で求めています。つまり... #ref(img56.png,nolink) とする事を繰り返します。 #ref(img57.png,nolink) これを繰り返し全ての非対角成分をゼロにし、対角化をするの... #ref(img58.png,nolink) この部分が行列Aを対角化する直交行列であり、列ベクトルは固... * Jacobi法の流れ [#vddf9a9e] まず、行列Aの非対角成分の中から絶対値が最大の成分を探しま... #ref(img59.png,nolink) のように、単位行列のpp成分をcosに、pq成分をsinに、 qp成分... #ref(img56.png,nolink) を考えたときに行列Aのpq、qp成分がゼロになるようにします。 #ref(img60.png,nolink) この4成分は展開すると上のようになり、pq、qp成分をゼロにす... #ref(img61.png,nolink) となり、 #ref(img62.png,nolink) #ref(img63.png,nolink) とおくと、 #ref(img64.png,nolink) また、 #ref(img65.png,nolink) とすると、 #ref(img66.png,nolink) #ref(img67.png,nolink) となります。 sign(αβ)とはαとβの積が正であればsinも正にな... #ref(img64.png,nolink) により、sinの符号がαβによって変化しています。 その他値が変わるのは次の式で表されるものだけで、それ以外... #ref(img68.png,nolink) * ソースコード [#ie265749] それではソースコードを見ていきます。ほとんど上の式を放り... for(i=0; i<eigenvectors->width; i++){ for(j=0; j<eigenvectors->height; j++){ if(i==j) eigenvectors->a[i*eigenvectors->width + j] ... else eigenvectors->a[i*eigenvectors->width + j] = 0; } } eigenvectorsを単位行列になるように初期化しています。これ... if(!(max = GetMaxvalue(matrix, &p, &q))) break; app = matrix->a[p*matrix->width + p]; apq = matrix->a[p*matrix->width + q]; aqq = matrix->a[q*matrix->width + q]; ここからループに入ります。 GetMaxvalueによって絶対値が最... alpha = (app - aqq)/2; beta = -apq; gamma = fabs(alpha)/sqrt(alpha*alpha + beta*beta); s = sqrt((1 - gamma)/2); c = sqrt((1 + gamma)/2); if(alpha*beta < 0) s = -s; 上の式の通り値を設定しています。 for(i=0; i<matrix->height; i++){ temp = c*matrix->a[p*matrix->width + i] - s*matrix->a[... matrix->a[q*matrix->width + i] = s*matrix->a[p*matrix->width + i] + c*matrix->a[q*mat... matrix->a[p*matrix->width + i] = temp; } for(i=0; i<matrix->height; i++){ matrix->a[i*matrix->width + p] = matrix->a[p*matrix->w... matrix->a[i*matrix->width + q] = matrix->a[q*matrix->w... } matrix->a[p*matrix->width + p] = c*c*app + s*s*aqq - 2*s... matrix->a[p*matrix->width + q] = s*c*(app-aqq) + (c*c - ... matrix->a[q*matrix->width + p] = s*c*(app-aqq) + (c*c - ... matrix->a[q*matrix->width + q] = s*s*app + c*c*aqq + 2*s... 式の通りなのですが、１つ目のfor文の1行目でわざわざtempに... for(i=0; i<matrix->height; i++){ temp = c*eigenvectors->a[i*eigenvectors->width + p] - s*eigenvectors->a[i*eigenvectors->width + q]; eigenvectors->a[i*eigenvectors->width + q] = s*eigenvectors->a[i*eigenvectors->width + p] + c*eigenvectors->a[i*eigenvectors->width + q]; ... eigenvectors->a[i*eigenvectors->width + p] = temp; } 直交行列Pの更新です。行列Aよりも処理が少ないのは、Pは右側... - &ref(main.c); - &ref(calculation.c); - &ref(calculation.h);

終了行:
RIGHT:寄稿：東條遼平 * 固有値と固有ベクトルを求める [#d575c2f8] Jacobi法(ヤコビ法)を用いて行列の固有値と固有ベクトルを求... * 固有値、固有ベクトルとは [#rdf1bce0] 任意のn次正方行列Aに対し、 #ref(img53.png,nolink) が成り立つときλを固有値、xを固有ベクトルという。 つまり、あるベクトルxに行列Aを掛けるとベクトルxが定数倍に... * Jacobi法について [#tb97b0dd] #ref(img54.png,nolink) となるとき、行列Aと行列Bの固有値は一致します。また、行列A... #ref(img55.png,nolink) と、対角化することができ、行列Λの対角成分は行列Λの固有値... このことから行列Aの固有値と固有ベクトルを求めるには、行列... Jacobi法ではこの直交行列Pを反復作業で求めています。つまり... #ref(img56.png,nolink) とする事を繰り返します。 #ref(img57.png,nolink) これを繰り返し全ての非対角成分をゼロにし、対角化をするの... #ref(img58.png,nolink) この部分が行列Aを対角化する直交行列であり、列ベクトルは固... * Jacobi法の流れ [#vddf9a9e] まず、行列Aの非対角成分の中から絶対値が最大の成分を探しま... #ref(img59.png,nolink) のように、単位行列のpp成分をcosに、pq成分をsinに、 qp成分... #ref(img56.png,nolink) を考えたときに行列Aのpq、qp成分がゼロになるようにします。 #ref(img60.png,nolink) この4成分は展開すると上のようになり、pq、qp成分をゼロにす... #ref(img61.png,nolink) となり、 #ref(img62.png,nolink) #ref(img63.png,nolink) とおくと、 #ref(img64.png,nolink) また、 #ref(img65.png,nolink) とすると、 #ref(img66.png,nolink) #ref(img67.png,nolink) となります。 sign(αβ)とはαとβの積が正であればsinも正にな... #ref(img64.png,nolink) により、sinの符号がαβによって変化しています。 その他値が変わるのは次の式で表されるものだけで、それ以外... #ref(img68.png,nolink) * ソースコード [#ie265749] それではソースコードを見ていきます。ほとんど上の式を放り... for(i=0; i<eigenvectors->width; i++){ for(j=0; j<eigenvectors->height; j++){ if(i==j) eigenvectors->a[i*eigenvectors->width + j] ... else eigenvectors->a[i*eigenvectors->width + j] = 0; } } eigenvectorsを単位行列になるように初期化しています。これ... if(!(max = GetMaxvalue(matrix, &p, &q))) break; app = matrix->a[p*matrix->width + p]; apq = matrix->a[p*matrix->width + q]; aqq = matrix->a[q*matrix->width + q]; ここからループに入ります。 GetMaxvalueによって絶対値が最... alpha = (app - aqq)/2; beta = -apq; gamma = fabs(alpha)/sqrt(alpha*alpha + beta*beta); s = sqrt((1 - gamma)/2); c = sqrt((1 + gamma)/2); if(alpha*beta < 0) s = -s; 上の式の通り値を設定しています。 for(i=0; i<matrix->height; i++){ temp = c*matrix->a[p*matrix->width + i] - s*matrix->a[... matrix->a[q*matrix->width + i] = s*matrix->a[p*matrix->width + i] + c*matrix->a[q*mat... matrix->a[p*matrix->width + i] = temp; } for(i=0; i<matrix->height; i++){ matrix->a[i*matrix->width + p] = matrix->a[p*matrix->w... matrix->a[i*matrix->width + q] = matrix->a[q*matrix->w... } matrix->a[p*matrix->width + p] = c*c*app + s*s*aqq - 2*s... matrix->a[p*matrix->width + q] = s*c*(app-aqq) + (c*c - ... matrix->a[q*matrix->width + p] = s*c*(app-aqq) + (c*c - ... matrix->a[q*matrix->width + q] = s*s*app + c*c*aqq + 2*s... 式の通りなのですが、１つ目のfor文の1行目でわざわざtempに... for(i=0; i<matrix->height; i++){ temp = c*eigenvectors->a[i*eigenvectors->width + p] - s*eigenvectors->a[i*eigenvectors->width + q]; eigenvectors->a[i*eigenvectors->width + q] = s*eigenvectors->a[i*eigenvectors->width + p] + c*eigenvectors->a[i*eigenvectors->width + q]; ... eigenvectors->a[i*eigenvectors->width + p] = temp; } 直交行列Pの更新です。行列Aよりも処理が少ないのは、Pは右側... - &ref(main.c); - &ref(calculation.c); - &ref(calculation.h);

ページ名: