概要
1、c++ では、行列は、多次元の配列でも一列に並んで格納されている。これをポインターを使用し、確認した。 更に、行列、配列を STL/array
と伴に、STL/vector も使用し、計算を行う予定である。
2、逆行列(余因子)
逆行列については、完了
・直接算出で、8行8列まで計算の実施(コード記載は6行6列まで、7行7列や8行8列は式が長いので未記載)
・連立一次方程式の解法( Gauss-Jordan 法)により、逆行列の算出。 (逆行列の本来の主目的に反して)
3、線形変換、固有値
線形変換は、当面完了とする。後日、追加して行く。
固有値は、対称行列の固有値算出法のJacobi法の計算のコード化を実施。
3行3列および4行4列の固有値を直接法にて、算出計算を実施。
5行5列、6行6列、7行7列の行列の要素を基にした固有多項式を作成した。
フレーム法の固有多項式の解をベヤストウ法で求める事を行っているが、ここには記載しない。
べき乗法での固有値算出を行っているが、ここには記載しない。
LR 分解での固有値と固有ベクトルの算出を行っているが、ここには記載しない。
Hessenberg化した行列を回転行列による QR 分解で、固有値算出のコードを記載した。
QR分解では、計算順序を考慮する事で、関数の低減により、繰り返し部の行数を低減させたコードにした。
programming 言語の c++ では、行列あるいは配列は、ポインターであると説明される。これを事例で確認した。 また、行列および配列が、どのように一列に配置されているかをポインターの順序で調査した
これの目的は、c++ の array は、最初に double a[2][3] のように、最初にデータ数の枠取りを行う必要がある。
このために、想定以上のデータが存在した場合には計算できない事になるので、データ数に実際上制限の無い
STL/vector で、行列・配列の計算が出来ないかのための準備でもある。
n 次行列の n 次行列式は、n! 個の項の和である。行列式値は、行と列を入れ替えても変わらない等々の性質を持つ。 n 次行列式 D の i
行と k 列を取り除いた n-1 次行列式を D(i,k) とする時に、 D の n-1 次小行列式と呼ぶ。(参考図書 1 ) この小行列式は、余因子行列と関係している。
上記 D の下線 ー は、D の上にあるべきであるが、書式の制約で表示できないので、下線の ー で示した。
今まで、参考図書は、各節の末尾に図書名のみを記載していたが、諸情勢より文中にも挙げる。
逆行列により連立一次方程式を解く事が出来る。また、行列式を利用し、連立一次方程式を解く方法に、Cramer の公式の手法がある。 詳細は、参考図書
1 に説明されている。なお、連立一次方程式の数値計算による解法として、 Gauss-Jordan 手法があり、本 Home page の第4章
数値計算 2 をベースに逆行列算出を下記 7、に示すのように行った。
行列 A = [ aik ] の逆行列 A-1 は、余因子行列 adj A と関係している。 A-1 = |A|-1 adj A となる。
また、AA-1 = A-1A= E である。逆行列の定義や余因子行列 adjA の展開式等々については、参考図書 1 をご
参照いただき、記載しない。行列式 A の値 |A| は det_A と表示している。 逆行列の算出を以下に示す。
1、余因子および逆行列の直接算出 コード (コード 2-1 )
計算コードは、配列の各要素ををそれぞれ一つの変数として、扱っている。
配列は、a[0][0] あるいは a00 のように、0 から始まるようにした。多くの参考図書のように a[1][1] あるい
は a11 ではない。
なお、計算結果を変数あるいは各余因子毎に出力するとなると、コードの行数が多くなるので、ポインタで
STL/vector に代入し、その vector により出力するようにした。例えば、ポインタ *ptr_a を定義し、
ptr_a=&aa[0][0] を上記の第一節に示すように、最初の参照位置としている。
コードでの各余因子行列は、参考図書の定義に従い、定式化を図り、それぞれの変数の積と和で算出した。
例えば、a11*a22 - a12*a21 である。a[1][1]*a[2][2] のように記載していない。定式化は気力を要する。
計算結果は、行列式 A と逆行列 adjA の積が、単位行列になる事を利用し、コードに示すように検証した。
1)、三行三列の行列式の逆行列
行列式の配列 a[3][3] の時の要素と変数を下記とする。
a00=a[0][0]、 a01=a[0][1]、 a02=a[0][2]、
a10=a[1][0]、 a11=a[1][1]、 a12=a[1][2]、
a20=a[2][0]、 a21=a[2][1]、 a22=a[2][2]。
この時の 余因子行列 hh[3][3] は、以下となる。ここで、hh[3][3] の [3][3] は 3 行 3 列である事を示す。
hh[0][0]= a11*a22-a12*a21、 hh[0][1]=-(a01*a22-a02*a21)、hh[0][2]= a01*a12-a02*a11、
hh[1][0]=-(a10*a22-a12*a20)、hh[1][1]= a00*a22-a02*a20、 hh[1][2]=-(a00*a12-a02*a10)、
hh[2][0]= a10*a21-a11*a20、 hh[2][1]=-(a00*a21-a01*a20)、hh[2][2]= a00*a11-a01*a10。
また、行列式 a の値 det_a は、次となる。
det_a=a00*a11*a22+a01*a12*a20+a02*a10*a21-a02*a11*a20-a00*a12*a21-a01*a10*a22;
これらにより、行列式 a の逆行列 a-1 は、余因子行列/det_a で、求める事が出来る。コード 2-1 。
なお、具体的な数値計算に用いた数値は、参考図書 1 の P.66 の演習問題 III の 1. i) に依った。
2)、四行四列の行列式の余因子行列と逆行列
行列式の配列 b[4][4] の時の要素と変数を下記とする。
b00=b[0][0]、b01=b[0][1]、b02=b[0][2]、b03=b[0][3]、
b10=b[1][0]、b11=b[1][1]、b12=b[1][2]、b13=b[1][3]、
b20=b[2][0]、b21=b[2][1]、b22=b[2][2]、b23=b[2][3]、
b30=b[3][0]、b31=b[3][1]、b32=b[3][2]、b33=b[3][3]。
この時の余因子行列 pp[4][4] は、以下となる。ここで、pp[4][4] の [4][4] は 4 行 4 列である事を示す。
この時の余因子行列の各要素の具体的な式については、コード 2-1 をご参照ください。
また、行列式 b の値 det_b についても、コード 2-1 をご参照ください。
これらにより、行列式 b の逆行列 b-1 は、余因子行列/det_b で、求める事が出来る。
なお、コードの具体的な数値計算に用いた数値は、参考図書 1 の P.67 の演習問題 III の 6. に依った。
3)、五行五列の行列式の余因子行列と逆行列
行列式の配列 c[5][5] の時の要素と変数を下記とする。
c00=c[0][0]、c01=c[0][1]、c02=c[0][2]、c03=c[0][3]、c04=c[0][4]、
c10=c[1][0]、c11=c[1][1]、c12=c[1][2]、c13=c[1][3]、c14=c[1][4]、
c20=c[2][0]、c21=c[2][1]、c22=c[2][2]、c23=c[2][3]、c24=c[2][4]、
c30=c[3][0]、c31=c[3][1]、c32=c[3][2]、c33=c[3][3]、c34=c[3][4]、
c40=c[4][0]、c41=c[4][1]、c42=c[4][2]、c43=c[4][3]、c44=c[4][4]。
この時の 余因子行列 qq[5][5] の各要素の具体的な式については、コード 2-1 をご参照ください。
ここで、qq[5][5] の [5][5] は 5 行 5 列である事を示す。
これらにより、行列式 c の逆行列 c-1 は、余因子行列/det_c で、求める事が出来る。
なお、コードの具体的な数値計算に用いた数値は、円周率 パイの数値列であり、逆行列の検証は、コード
に示す検証計算で、E になる事を確認した。この検証のコードは、3次、4次の場合と同じであり、3次、4次
の参考の計算例でも、E となるので、誤りはない。
4)、六行六列の行列式の余因子行列と逆行列
この行列式はコードで、コード 2-1 に追加した。一式が長い。
5)、七行七列、八行八列の余因子行列と逆行列
余因子行列と逆行列の各成分算出する上記のような式の計算コードを作成し、計算検証を行った。しかし、
各成分の式が非常に長く、数も 36、49 と多いので、コードをここに記載しない。
2、連立一次方程式の解法( Gauss-Jordan 法)による逆行列 コード(コード 2-2 )
この Gauss-Jordan 法で、逆行列の各成分あるいは要素を算出した。更に、行列式値も算出した。
計算に用いたコードは、c++ 数値計算 2 第四章 連立一次方程式の 1、Gauss-Jordan で用いたコードを
ベースにしている。連立一次方程式の解法で、逆行列を算出する手法は、目的と手段が入れ替わったと考え
られが、STL/vector を用いた事で汎用性が上がった。
計算の概要
・入力 全ての正則行列の逆行列算出を STL/vector で行うために、行列式 a の成分 a00、a01、a02、、
、a10、a11、a12、、、a20、a21、a22、、、を行列の形式での入力としないで、各成文を、
上記の順で一行で入力するようにしている。
この結果、行列式の各成分が、一列となる csvファイルで、コードでは、" data_0.csv " が入力
ファイルになる。これで、STL/array のように、予め、領域の枠取りを行う作業が無くなり、
3行3列、4行4列、、、に対して、それぞれコードを準備する必要がない。
読み込んだ入力データは、" in_det.csv " で確認できる。
・出力 行列式値および余因子行列と逆行列 B の算出結果が、各成分 b00、b01、b02、、、b10、b11、
b12、、、b20、b21、b22、、、の順に 各成分が一列に、 " inverse_determinant.csv
"
ファイルに出力される。
・計算結果の検証 行列と逆行列(この時は連立一次方程式の解の行で、転置行列の列ではない)の積の
対角要素が 1.0 になる事を "result_check.csv " ファイルで確認出来る。
また、参考計算の行列式の要素値は、上記の六行六列と同じで、同じ結果を得ている。
・行列値の算出に、連立一次方程式の行列式の要素を示す vector の s を用い、計算を簡素にした。
上記の 1、から 6、( 3 行 3 列から 8 行 8 列)までの定義にしたがった計算値と合っている事を確認。
なお、行列式値算出に 4 つの手法の計算を示した。
参考図書
参考図書 1、古谷 茂著:新数学シリーズ 行列と行列式(昭和45年3月 増補第29刷) 培風館
なお、スミルノフ高等数学教程も参考にしている。両者とも、線形変換が記載されている。
感想 コードは、c++の特長を生かした数値計算にはまだなっていない。更なる工夫が必要である。
また、文中にコードそのもの式の記載を行っていたが、泥臭すぎたので、式を削除した。
copyright©2009-20xx Nakadaka&Kobaru all rights reserved.
当 Page から他の page へ
Index c++から風と翼の計算へ
0、c++ 文字(半角/全角)、フォルダ生成、csv読み込み・書き込み、データ処理
1、c++ STL適用例
2、c++11調査、番外(多計算並列)
3、c++ 数値計算 1 三次・四次方程式、方程式求根、微分方程式、 最小二乗法、フーリエ変換
4、c++ 数値計算 2 連立一次方程式、補間法: Spline、Lagrange、Aitken、直交関数法
5、c++ 数値計算 3 逆ラプラス変換、移動則、Bode Plot
6、c++ 数値計算 4 Feedback回路、Euler変換(62次まで)、二項係数
7、c++ 数値計算 5 直交関数: Legendre(34次)、Chebyshev(35次)、Bessel 多項式級数展開:Legendre、Chebyshev
8、c++ 数値計算 6 数値積分(連続系): Newton-Cotes、Legendre、Gauss-Chebyshev
9、c++ 数値計算 7 数値積分(離散系): 不等間隔データ(Newton-Cotes法)
10、c++ 数値計算 8 行列式(余因子、 逆行列 )、行列(線形変換、固有値)
11、c++ 数値計算 9 積分方程式
12、c++ 数値計算 予備
13、c++ 技術計算 1 珪素鋼板 BH特性、BH特性の補間
14、c++ 技術計算 2 曲線の交点算出(Function文とClass)
15、予備
16、予備
17、よもやま話と時代の位相
18、予備
19、予備