DelFEMで行列の固有値を求める
DelFEMで行列の固有値を求めてみます。
DelFEMのサンプルを見ると連立方程式を解くdeterministic problem(決定問題？）は扱っていますが、固有モード解析の例が見当たりません。
しかしながらライブラリの中を探してみるとMinimumEigenValueVector_InvPower (delfem/eigen?lanczos.h) というのがありました。
簡単な例で試してみます。
※ソースコードは自作のC#用ラッパーライブラリDelFEM4Net(URL: http://code.google.com/p/delfem4net/ )を使用しています。
A = [ 1 1/2]
[1/2 1]
の固有値を求めます。
この行列の固有値は、3/2と1/2です。
using DelFEM4NetCad;
using DelFEM4NetCom;
using DelFEM4NetCad.View;
using DelFEM4NetCom.View;
using DelFEM4NetFem;
using DelFEM4NetFem.Field;
using DelFEM4NetFem.Field.View;
using DelFEM4NetFem.Eqn;
using DelFEM4NetFem.Ls;
using DelFEM4NetMsh;
using DelFEM4NetMsh.View;
using DelFEM4NetMatVec;
using DelFEM4NetLsSol;
Console.WriteLine("CEigenLanczos.MinimumEigenValueVector_InvPower");
double[,] srcMat = new double[2, 2]
{
{1, 0.5},
{0.5, 1}
};
// DelFEMの固有値解析ルーチンを使用
// リニア方程式ソルバ
CLinearSystem ls = new CLinearSystem();
// プリコンディショナ―
CPreconditioner_ILU prec = new CPreconditioner_ILU();
// 行列の０でない要素のある場所を指定し、インデックスを格納する
uint nblk = (uint)srcMat.GetLength(0);
ls.AddLinSysSeg(nblk, 1);
// 行列の要素のインデックスを格納するインスタンスを生成
CIndexedArray crs = new CIndexedArray();
crs.InitializeSize(nblk);
using (UIntVectorIndexer index = crs.index)
using (UIntVectorIndexer ary = crs.array)
{
index[0] = 0;
// 行を走査
for (int col = 0; col < nblk; col++)
{
// 列を走査
for (int row = 0; row < nblk; row++)
{
// ここではn x n の行列の要素すべてを格納できるようにします。
// 本当は非０の場合のみ追加するものだと思われます。
// 対角要素は除きます。
if (col != row)
{
ary.Add ( (uint)row );
}
}
index[col + 1] = (uint)ary.Count;
}
}
// 作成したインデックスをリニアシステムの対角行列？に登録
ls.AddMat_Dia(0, crs);
//prec.SetFillInLevel(1);
// ILU(0)のパターン初期化
prec.SetLinearSystem(ls);
// 行列の要素を格納する場所を取得
CMatDia_BlkCrs_Ptr femMat = ls.GetMatrixPtr(0);
// 行列へのマージ開始
ls.InitializeMarge();
//uint nblk = femMat.NBlkMatCol();
// 行列を要素を追加します。
System.Diagnostics.Debug.Assert(nblk == femMat.NBlkMatCol());
Console.WriteLine("srcMat merge");
for (uint iblk = 0; iblk < nblk; iblk++)
{
DoubleArrayIndexer ptr = femMat.GetPtrValDia(iblk);
System.Diagnostics.Debug.Assert(ptr.Count == 1);
ptr[0] = srcMat[iblk, iblk];
Console.WriteLine(" ( " + iblk + " ) = " + ptr[0]);
uint npsup = 0;
ConstUIntArrayIndexer ptrInd = femMat.GetPtrIndPSuP(iblk, out npsup);
DoubleArrayIndexer ptrVal = femMat.GetPtrValPSuP(iblk, out npsup);
System.Diagnostics.Debug.Assert(ptrInd.Count == ptrVal.Count);
for (int i = 0; i < ptrVal.Count; i++)
{
ptrVal[i] = srcMat[iblk, ptrInd[i]];
Console.WriteLine(" ( " + iblk + ", " + ptrInd[i] + " ) = " + ptrVal[i]);
}
}
// 更新ベクトル?
CVector_Blk_Ptr updatePtr = ls.GetUpdatePtr(0);
System.Diagnostics.Debug.Assert(nblk == updatePtr.NBlk());
System.Diagnostics.Debug.Assert(1 == updatePtr.Len());
for (uint iblk = 0; iblk < updatePtr.NBlk(); iblk++)
{
double val = 1.0;
updatePtr.SetValue(iblk, 0, val);
//System.Console.WriteLine("update( " + iblk + " ) = " + updatePtr.GetValue(iblk, 0));
}
// マージ終了
double res = ls.FinalizeMarge();
//Console.WriteLine("Residual : " + res );
// プリコンディショナ―にリニアシステムの(行列の？)値を設定してILU分解を行う
bool ret = prec.SetValue(ls);
//Console.WriteLine("prec.SetValue(ls) ret = " + ret);
// 行列の固有値を取得する
// 最大値が取得できるようです。
double conv_res_lssol = 1.0e-6;// ICCG法の収束基準の相対残差
uint itr_invp = 2000;// 逆べき乗法の最大反復回数
uint itr_lssol = 2000;// ICCG法の最大反復回数
int iflag_conv = 0; // 0:正常に終了 1:ICCGが収束しなかった
IList aIdVec = new List(); //何をセットすればいいのか?
/*
uint ntmp = ls.GetTmpVectorArySize();
for (uint itmp = 0; itmp < ntmp; itmp++)
{
aIdVec.Add(itmp);
}
*/
double eigenVal = CEigenLanczos.MinimumEigenValueVector_InvPower(
ls,
prec,
aIdVec,
itr_invp,
itr_lssol,
conv_res_lssol,
out iflag_conv);
// iflag_conv (0:正常に終了 1:ICCGが収束しなかった)
Console.WriteLine("////// Result of CEigenLanczos.MinimumEigenValueVector_InvPower //////");
Console.WriteLine("eigenVal = " + eigenVal);
Console.WriteLine("iflag_conv = " + iflag_conv);
for (uint iblk = 0; iblk < updatePtr.NBlk(); iblk++)
{
System.Console.WriteLine("eigenVec( " + iblk + " ) = " + updatePtr.GetValue(iblk, 0));
}
Console.WriteLine("press any key");
Console.In.Read();
実行結果
メソッドの名称は最小固有値ですが、最大値が求まりました。
InvPowerとあるのでこれでいいのだと思います。
これはエラーケースでした。逆べき乗法（逆反復法？）では固有値の最小値が求まるはずです。試しに行列の対角要素を少しずらすとうまく最小値を取得できました。
// 固有値の最大値が求まる → これはエラーケースのようです。
/*
double[,] srcMat = new double[2, 2]
{
{1.0, 0.5},
{0.5, 1.0}
};
*/
// 対角の値を少しずらすと最小値が求まります。
double[,] srcMat = new double[2, 2]
{
{1.0, 0.5},
{0.5, 0.999999999999999}
};
なお、固有ベクトルは
{1, 1}T (λ = 1.5に対して)
{1, -1} T (λ= 0.5に対して)
(λは固有値)
なので、こちらも問題なく取得できています。
