MATLABで学ぶ固有値問題の解析法　べき乗法（行列反復法）

2022.02.23

はじめに

MATLABで固有値問題を解くときは「eig」や「eigs」を使って求めていると思います。これらの関数で簡単に固有値問題が解けてしまうので、最近は「固有値問題の解法」についての知識がある人が少なくなっています。

モード解析やFEMでの振動解析を生業としている方には「固有値問題の解法」についての知識があったほうがいい、というのが私の考えです。

もし、「CAEとしてFEMをしています」という人（CAEのプロ）が有限要素法の内挿関数、グラスアワーモード、k6rotなどを知らない場合、そのCAEの解析結果や分析結果を心底信用できますか？

私は信用できないんですよね。。。
（誰だって間違いがありますが、間違えの可能性が（知識のある人よりも）少し高いのではないかと思いつつ、解析結果をレビューしている気がします。）

解析ソフトを使えば、高校生でも解析はできます。
ただ、その解析結果が妥当なのか、確認したい物理現象を考察するのに十分か、などについては周辺の深い知識が必要となる場合があります。

「固有値問題の解法」についても同様で、ある程度知識を持っていてほしいと考えています。そこで、本記事ではべき乗法（行列反復法）という固有値問題の解法を紹介します。

 

べき乗法（行列反復法）

ざっくり理論を紹介します。おなじみの式(1.1)を式(1.2)のように変形します。

$$
[K](x)=ω^2[M](x) \tag{1.1}
$$

$$
[A_1](x)=λ(x) \tag{1.2}
$$

べき乗法はこの固有値問題の最高次の固有値を求める方法です。
任意の初期値ベクトル(xo)を繰り返し計算により、真の解(x)に近づける方法です。

数式を用いた解説は大学の講義PDFを参考に頂けると幸いです。

https://www.cspp.cc.u-tokyo.ac.jp/hanawa/class/spc2016s/sp20160607-1.pdf

http://www.math.ritsumei.ac.jp/yasutomi/jugyo/Numerical_Analysis/2009/note15.pdf

MATLABなどのコードを読める方は、文章で説明するよりもプログラムを見たほうがわかりやすいと思いますので、下記の例題で解説いたします。

MATLABを用いた例題

下記のような質量行列[M]と剛性行列[K]があったとします。

clear all

[K]=eval_Kmatrix(ones(1,4)*1000);
[M]=eval_Mmatrix(ones(1,4));

A=M\K;

まずは、MATLABのeigsで固有値解析してみましょう。

[V2,D2] = eigs(A,ndof,'lm')
D2=diag(D2)

すると、下記のような答えになりますね。

V2 =

   -0.4285    0.6565    0.5774   -0.2280
    0.6565   -0.2280    0.5774   -0.4285
   -0.5774   -0.5774   -0.0000   -0.5774
    0.2280    0.4285   -0.5774   -0.6565


D2 =

   1.0e+03 *

    3.5321
    2.3473
    1.0000
    0.1206

 

さて、べき乗法での解法は下記です。

clear all

[K]=eval_Kmatrix(ones(1,4)*1000);
[M]=eval_Mmatrix(ones(1,4));
ndof=4;
A=M\K;
X0=[1 2 3 4].';
lambda=zeros(ndof,1);
V=zeros(size(A,1),ndof);




for ii1=1:ndof
    eps=10^-5;
    dm=1;
    X=X0;
    while abs(dm) >eps

        xs=sqrt(sum(X.^2));
        X=X/xs;

        X1=A*X;

        rm=sum(X.*X1);
        xs=sqrt(sum(X1.*X1));
        X1=X1/xs;

        for ii2=1:ndof
            if X(ii2)==0
                dm=X1(ii2)-X(ii2);
            else
                dm=(X1(ii2)-X(ii2))/X(ii2);
            end
        end
        X=X1;

    end


        lambda(ii1)=rm;
        V(:,ii1)=X;

        AA=X*X';
        A=A-rm.*AA;
end
V
lambda

結果は下記です。

V =

   -0.4285    0.6565   -0.5774    0.2280
    0.6565   -0.2280   -0.5773    0.4285
   -0.5774   -0.5773    0.0000    0.5774
    0.2280    0.4285    0.5773    0.6565


lambda =

   1.0e+03 *

    3.5321
    2.3473
    1.0000
    0.1206

どちらも同じですね。

functionファイル

なお、このとき文末の「functionファイル」をカレントディレクトリに保存しておいてください。

function [K]=eval_Kmatrix(k_vec)

K=zeros(length(k_vec));
for ii1=1:length(k_vec)
if ii1==1
K(ii1,ii1)=k_vec(ii1);
else
K(ii1-1:ii1,ii1-1:ii1)=K(ii1-1:ii1,ii1-1:ii1)+[1 -1;-1 1]*k_vec(ii1);
end
end

end

function [M]=eval_Mmatrix(m_vec)

M=diag(m_vec);

end

 

 

コメント