非対称特異行列のクリティカルな左固有ベクトルの力学的意味付けとその例説

日本機械学会論文集（A 編） 原著論文 No.2013-JAR-0288

非対称特異行列のクリティカルな左固有ベクトルの

力学的意味付けとその例説

＊

Fumio FUJII

*1

, Yuki YAMAKAWA, Yoshihiro INOUE,

Yasuko MIHARA and Takaya KOBAYASHI

In mathematical problems and mechanical engineering, there are a number of examples, in which a non-symmetric Jacobian matrix is involved in solution of simultaneous linear equations. The left and right eigenvectors of such non-symmetric matrices are in general complex and must be well discerned to each other. Their properties and practical meaning are, however, hardly discussed in engineering applications. Especially when the non-symmetric matrix is singular, the critical left eigenvector corresponding to null eigenvalues is of increased significance in the examination of the solvability of the problem. The present paper describes and interprets the substantial role of the critical left eigenvector of the non-symmetric singular matrix in mechanics. Model examples in applied mathematics, solids, structures and rigid bodies will illustrate the meaning of the critical left eigenvector, when the singularity is unavoidable in the problem to be solved. The discussion will be then extended to the critical left singular vector of a rectangular matrix.

Key Words : Computational Mechanics, Structural Analysis, Numerical Analysis, Finite Element Method, Nonlinear

Problem 1. 緒 言 固体，構造，流体，制御，ロボット，電気回路などの工学分野の非線形応答の解析は，次のような線形化され た方程式の連続的な更新と，その解法に帰着されることがある(1)~(6)_： = Jx b （1） ここで一般に， T (≠ ) J J ：N×Nの実非対称行列， x ：解ベクトル，b：右辺項である．また本稿では，Jが例 えば分岐座屈点の前後で更新の結果，正則行列から特異行列となる場合を想定している． detJ=0 （2）

藤井

文夫

*1

，山川

優樹

*2

，井上

吉弘

*1

三原 康子

*3

，小林 卓哉

*4

Interpreting the Critical Left Eigenvector of a Non-Symmetric Singular Matrix in Mechanics

*1 _{Gifu Univ., Dept. of Mechanical Engineering,} 1-1 Yanagido, Gifu City 501-1193, JAPAN

* _{原稿受付 2013 年 4 月 2 日} *1 _{正員，岐阜大学工学部機械工学科（〒501-1193 岐阜県岐阜市柳戸 1-1）} *2 _{東北大学大学院工学研究科土木工学専攻（〒980-8579 宮城県仙台市青葉区荒巻字青葉 6-6-06）} *3 _{（株）メカニカルデザイン（〒182-0024 東京都調布市布田 1-40-2）} *4 _{正員，（株）メカニカルデザイン} E-mail: [email protected] 79 巻 808 号 （2013-12）

このような更新に伴い特異点でゼロとなる固有値λ_cをクリティカルな固有値，対応する左右の固有ベクトル c c (φ θ, )をクリティカルな固有ベクトルと呼ぶことにする．さて線形化された方程式（1）の求解においては，特 異性のため解けないとするのではなく，特異行列の数理構造を考察し，右辺項bを操作することで，求解にむけ ての状況を積極的に作り出すことが，新技術の開発ヒントとなることがある． 著者らはすでに，数万から数十万自由度の規模の非線形有限要素解析において登場する対称剛性行列 （_K=KT_{）について，固有値解法を経由することなくLDL 分解法の副産物として，特異性に絡むクリティカ}T ルな固有ベクトル（座屈モード）を抽出できる実用的な計算手法を提案した(7)~(10)_{．対称行列においては，左固有} ベクトルと右固有ベクトルとの間の区別はない．しかし非対称行列（ _≠ T J J ）においては，左固有ベクトルφjと 右固有ベクトルθj（j= 1, ,N）は一般に，複素ベクトルとなり，互いに異なるモードとなる．文献(11)(12)におい て著者らは，クリティカルな左右の固有ベクトルもLDU分解法の副産物として抽出できるアルゴリズムを提案 し，非線形有限要素法による材料分岐解析への適用も行った．しかし左右の固有ベクトルのもつ特質や意味につ いて，工学的応用において着目されることは稀である．そこで本稿では，特にクリティカルな左固有ベクトルの 工学的な意味付けを考え，具体的な力学モデルを使って例説する． 2. 左右固有ベクトルと可解条件 線形代数学を簡単にレビューする(13)~(19)_{．まず非対称行列} J（N×N）は 固有値分解が可能であるとして，左 固有ベクトルφjと右固有ベクトルθjを，次のように定義しておく． T T j J =λj j φ φ （3） および j =λ j j Jθ θ （4） である．式（3）および式（4）において，j= 1, ,N であり，λjは左右の固有ベクトルが共有する固有値であ る．さらに式（3）と式（4）の両辺を転置すると，それぞれ次のようになる． T j =λ j j J φ φ （5） および T T T jJ =λj j θ θ （6） これよりJの左固有ベクトルφjはJTの右固有ベクトルであり，Jの右固有ベクトルθjはJTの左固有ベクト ルでもある．したがって，Jの左右の固有ベクトルの工学における現実的な意味は，Jと T J との間で力学的双 対性が見出せる場合にはより鮮明となる．この点にも本稿では着目する． 非対称行列の固有解(λ φ θj, j, j)は一般に，複素数・複素ベクトルとなる．しかしλc= の左右の固有ベクト0 ル(φ θc, c)は，実ベクトルとして選ぶことができる(12). 異なる固有解の組み合わせでは，左右の固有ベクトルに ついては，つぎのような直交性があり， T 0 for j k = j≠k φ θ （7） 同じ固有解の組では，つぎのように規格化しておく．

T 1 for 1, , j j = j= N φ θ  （8） ただしベクトル成分の比を強調するために，固有ベクトルを必ずしも規格化せずに提示することもある．また Jについては，次のようなスペクトル展開が可能である． T 1 N j j j j λ = =

∑

J θ φ （9） ここで， T j j θ φ は射影行列であり，式（8）より T T T j j j j = j j θ φ θ φ θ φ （10） が成立する．またJが正則行列である場合に限り形式的に， T 1 1 N j j j j λ − = =

∑

J θ φ （11） が成立する．そしてこのとき，式（1）の解は一意的に，次のように決まる． T 1 N j j j j= λ =

∑

b x φ θ （12） すなわちJが正則行列のとき，右固有ベクトルθ jは，解を構成する線形独立なベクトルとなる．θ jの重みを 決定する分数表現では，分子で左固有ベクトルφjと右辺項bとの内積が，分母では固有値λjが関与してくる．J がλc→ となるような特異行列のとき，クリティカルな左右の固有ベクトルを0 (φ θc, c)として，解が存在するた めの条件は，次のようになる(15)_． T c b=0 φ （13） これより可解条件として，クリティカルな左固有ベクトルであるφcと右辺項bとの間の直交条件が必要となり， c θ の重みは，0 0の不定形となる．この不定値を決めるためには解くべき問題に応じて一般に，追加的な条件式 や支配方程式の高階導関数が必要となる(1)_{．ここまでが数学的な意味付けである．} 3．簡単な例題における左固有ベクトルの役割の例説 本稿の主旨を，まず簡単な平面座標系の変換の例題で示しておく．座標の( , )p q ←( , )s t なる順方向の簡単な 1 次変換を考える． 1 3 2 6 p s q t   ₌               （14）

変換行列の第 1 行の 2 倍が第 2 行であるため，q=2s+6t=2(s+3 )t =2pより，すべての写像点( , )p q は直線 2 q= pの上に点在する．図 1 に代表的なサンプル点( , )s t に対する写像点( , )p q を示す．変換行列が特異行列であ っても順方向の変換は成立し，変換後の点を実際にプロットすることで点列が並ぶ直線を特定することができる． ちなみに式（14）の中にある特異行列の固有値は，λ1=λc= および0 λ2 = であり，対応する左右の固有ベク7 トルは T c= −( 2,1) φ , θc= −( 3,1)Tおよびφ2 = − −( 1, 3)T, θ2 = − −( 1, 2)Tである． しかしこの直線の方程式q=2pは，逆方向の変換( , )p q →( , )s t ではどのような操作から抽出することができる のであろうか．変換行列は特異行列であるため一意的な逆変換は存在しないが少なくとも形式的には，与えられ た右辺項( , )p q に対して解( , )s t を求める次のような連立方程式の解法に帰着させることはできる． 1 3 2 6 s p t q      ₌             （15） この場合には，すべての右辺項( , )p q に対して解( , )s t が存在するわけではなく，左固有ベクトルである T c= −( 2,1) φ と右辺項( , )p q との直交条件から， 2p− = の場合に限り不定解q 0 ( , )s t が存在する．すなわち特異 行列のクリティカルな左固有ベクトルに着目することにより，式（15）が解をもつような右辺項の設定操作によ り，特異行列であるにも関わらず逆変換( , )p q →( , )s t が可能となる状況を生みだすことができる．これがクリ ティカルな左固有ベクトルの意味付けとなる．したがって式（13）で示したように，問題として成立するか否を 判定する左固有ベクトルφcの持つ役割は本質的に重要となる．ちなみに右固有ベクトルθc= −( 3,1)Tは式（15） において，右辺項( , )p q =(0, 0)に対する同次解( , )s t を形成する． 4. 反対称行列のクリティカルな左固有ベクトル この章では反対称行列について例説する．章全体の展望を表 1 にまとめておく．表 1 を掲載する理由は，異な る例題でも数理的に同じ本質を共有するベクトル間の対応を明示し，ひいては各例題における式（13）の可解条 件を整理しまとめて展望するためである． ( , )p q

⇐

( , )s t

Fig.1 Linear coordinate transformation

Table 1 Skew-symmetric matrix and cross-product matrices in Chap. 4 summarized

Section # Problem Setup Matrix Critical Left Eigenvector To be solved Specified Solvability or Orthogonality 4.1 Linear Algebra T = − J J J Skew- Symmetric Matrix c φ Jx=b J, b φc, b 4.2 Cross Product c= ×a b A Cross- Product Matrix a _Ab=_{c A, c} a , c 4.3 Moment m= ×r f R r Rf = m _{R, m} r , m 4.4 Velocity v= ×ω r W ω Wr=v W , v ω, v

4・1 反対称行列の数理 三次元剛体力学や連続体力学で登場する反対称行列は，左固有ベクトルの意味付けを議論する観点からは好例 であり，その特質をまとめておく．反対称行列Jの定義は，次のようである． T = − J J （16） またはJの対称部分が消滅することから T + = J J O （17） でもある．一般に反対称行列のランクは偶数で，固有値はゼロまたは純虚数である(19)_{．ここで式（16）から} T T

detJ =det (−J )= −( 1)NdetJ = −( 1)NdetJ （18）

となる．ここで T

detJ =detJを活用した．したがってNが奇数のとき，式（18）より detJ=0となり，その左

右のクリティカルな固有ベクトル(φ θc, c)については， T T T T T c J= → − c J = → J c = φ 0 φ 0 φ 0 （19） および T T T c= → − c= → c = Jθ 0 J θ 0 θ J 0 （20） である．したがって奇数次の反対称行列では，φcおよびθcの区別は消滅し，両者のモードは一致する．特にN=3 の反対称行列は連続体力学でも高い頻度で登場し，次節でも示すように三次元空間でのベクトル積を記述できる ため（図 2 と図 3），クリティカルな左固有ベクトルの意味づけを議論する題材としても重要である．

(a) Moment m (b) Angular velocity v Fig.2 Vector product Fig.3 3D rigid body mechanics

4・2 ベクトル積の外積行列による表現

ベクトル c を，ベクトル a およびbの外積として定義する（図 2）：

= ×

a とbが極性ベクトルのとき， c は軸性ベクトルとなる．式（21）の右辺は，bの左側から作用するa× •( )の演 算を外積行列Aを使って，次のように書き換えることができる． = c Ab （22） ここに外積行列として， 3 2 3 1 2 1 0 0 0 a a a a a a −     =_ − _ −    A （23） を導入し，また

(

)

T 1 2 3 a a a = a （24） である．外積行列Aは反対称行列であり， a をAの双対ベクトルと呼ぶ(4)．そして， T T T (− ) − − → = and = A a = A a = Aa = a × a =0 a A 0 Aa 0 （25） であることから，双対ベクトル a はAのクリティカルな左右の固有ベクトルでもある． 幾何学的な別解釈として式（22）を，a に対して有限の回転操作を与えて，bに変化させる回転ベクトルが c で あるとして解釈できる（図 2）．そこで式（22）を書き換えて， a に c なる回転操作を与えて生成されるベクトル bを求める問題を，次のような求解問題に帰着させて考えてみる． = Ab c （26） このとき detA=0であるので，回転操作の結果としての b が存在するためには，既定すべき c には拘束が課せら れ，Aのクリティカルな左固有ベクトル a と c とが直交する場合に限り，回転操作の結果としてのbが可能とな る． T 0 = a c （27） これがクリティカルな左固有ベクトル a の工学的な存在意義である．すなわち三次元ベクトル a に回転ベクト ル c を作用させてbを生成する有限回転運動を式（26）で記述する限り，bの求解問題としての式（26）が成立 するためには式（27）が不可欠となる．bに課せられた条件は特にないが，結果的には a とbが，回転ベクトル c を法線ベクトとする面内にとどまることになる．ところで， × = − × a b b a （28） でもあるから， = − Ab Ba （29）

も成立する．Bはbの外積行列である．本節では外積行列の双対ベクトル，すなわちクリティカルな左右の固有 ベクトルの意味づけを考えた．ベクトル積は次節でも示すように，剛体の回転平衡や回転運動を議論する際にも 使われる． 4・3 位置ベクトルと力ベクトルの外積としてのモーメント（図 3 (a)） 古典力学では，点 O から伸びる位置ベクトルrがあり，その先端 P に作用する力のベクトル f が生み出すモー メント m は，一般に次のように定義される． = × m r f （30） ベクトル積として生み出される軸性ベクトル m は結果的に，rと f の両方に直交する．しかし逆に，アーム長 r（または外積行列R）に対して，結果としての m を生み出すような力 f を求解する式は，次のようである． = Rf m （31） これより力学的に可能な解 f の存在条件は，rと m とが直交することである．すなわちrを法線ベクトルと する平面内に m があるような条件設定が m に課せられる．この条件設定においても，Rのクリティカルな左固 有ベクトルであるrが役割をもつ．なお節点荷重による有限要素の平衡を考えるときのように，点Oと点 P の間 に充填された材料の分布形状は，位置ベクトルrのような直線である必要はなく，図 3 (a) の点線のように両端 は任意の形状のアームで連結されていても構わない．両端の座標のみが重要である． 4・4 位置ベクトルと回転速度ベクトルの外積としての速度ベクトル（図 3 (b)） 前節の力学式と酷似した式に，運動学の式がある．線分rが瞬間的にωなる角速度運動を受けたとき，rの先 端の速度 v は，次のようである． = × v ω r （32） この外積計算ではωが軸性ベクトルで， r は極性ベクトルなので，速度ベクトル v は極性ベクトルとなる．前 節の例題との間では( , , )m r f ⇔( , , )vω r の対応があり，この対応を表 1 と図 3 でも強調してある．ω（または外 積行列W ）と v を与えて，未知量 r を求める問題は，次のような方程式の解法に帰着される． = Wr v （33） ここでも，ωと v を任意に規定することはできず，ωを法線ベクトルにもつ面内（図 3 (b) の円内）に v が収 まる必要がある．すなわち角速度ωに対して速度 v が直交するように問題設定をしない限り，解 r が存在する有 意な工学問題として成立しない．これがW のクリティカルな左固有ベクトルωのもつ意義である．なお解 r はω を法線とする面内にある必要はない（ωと r が定義する平面の法線が v となる）． 5．アームの特異姿勢 JとJTが双対問題を記述する場合には，双対問題が切り替わる度に左右の固有ベクトルの役割が交換される ため，より明確な工学的解釈が可能となる．例として図 4 にある単位長さの剛体棒（L=1.0）を 2 本直列に連結 したロボットアームを考える．関節点 A, B における折れ角( , )α β を用いて先端点 C の座標( , )x y を求めると，次 のようになる．

cos cos( ) sin sin( ) x y α α β α α β + +   ₌     _{+ +}      （34） さらに時間微分にドット（・_{）を導入し，式（34）の両辺を速度表現し左右辺を交換する．するとアーム先端の} 並進速度( , )x y  を可能とする節点角速度( , )α β を求める求解問題は，次のようになる．

Fig.4 Kinematics and equilibrium of a link arm

x y αβ       α_β =     J _   （35） ここで

sin sin( ) sin( )

cos cos( ) cos( )

αβ = − _αα− _{α β}α β+ − _{α β}α β+ _ + + +   J （36） である．このヤコビアン行列の行列式はdetJ=sinβ で，β =0のとき特異行列となる．そこで具体的な特異姿 勢としてtanα =0.75, α =36.87, β =0を仮定すると， 1.2 0.6 1.6 0.8 − −   =     J （37） と なる ．こ のJ の 固有 値はλ1= −0.4お よびλ2=λc = であり，左右の固有ベクトルは0 φ1= − −( 2, 1)T , T 1= −( 3, 4) θ および T c= − −( 4, 3) φ , θc=(1, 2)− Tである．

(a) Velocity vector of arm tip C (b) Angular velocity of joints A & B Fig.5 Meaning of the critical left and right eigenvectors of J in kinematics

αβ = J Jの時の式（35）において関節角速度( , )α β を得るための設定条件として，クリティカルな左固有ベクト ルφcを参照して，右辺項に対して次式が課せられる．

(

4 3

)

x 0 y   − −  =     （38） これより，アーム先端速度はアーム軸線方向と直交する場合に限り，解としての関節点角速度を得ることができ る（図 5 (a)）．これ以外のアーム先端速度に対しては，解として関節の角速度は存在しないので，工学問題とし ては成立しない．このようにアーム先端の速度を見出す際に，左固有ベクトルφcは運動学的に重要な役割をなす． 一方，クリティカルな右固有ベクトルθcをJの右からあてがうと，次のようになる． 1.2 0.6 1 0 1.6 0.8 2 0 − −       =      _{ −}       （39） すなわち，関節に右固有ベクトルのモードが示す回転速度を与えてアームの運動を制御すると，特異姿勢にある アーム先端は瞬間的に静止を保ったままの状態となる（図 5 (b)）．これが同次解でもある右固有ベクトルθcの運 動学的な解釈である． 次に問題を裏返して，同じ特異姿勢にあるアームの平衡問題を考える（図 6）．関節点に駆動力(mA,mB)を与 えた場合に先端で許容できる荷重( , )p q を求める問題は，次のような式の解法に帰着される． A T B p m q m     =         J （40）

(a) Joint torque distribution (mA,mB) (b) Tip loading for null torque Fig.6 Meaning of the critical left and right eigenvectors of JT in equilibrium

このときJの左右の固有ベクトルが逆転する結果， T J のクリティカルな左右の固有ベクトルはそれぞれ，θcお よびφcとなる. 式（40）の解法において現実的な解( , )p q が求まるためには，今度はθcを参照して，次のような 駆動力を与えなければならない．

(

)

A B 1 2 m 0 m   −  =   （41） 式（41）を満たす駆動力の比を，節点角速度に例えて幾何学的に図示したのが図 6 (a) である．式（41）を満た さない駆動力は，先端荷重に対して平衡を保つことはできないことになるが，これは真直ぐに伸びたアームの特 異姿勢において，先端から各関節までのアーム長の比を考慮すると力学的に当然である．同じことの数理表現が c θ を用いた式（41）となる． T J の右固有ベクトルφcは式（40）の同次解であるので，JTの右からφcをあて がうと右辺項はゼロベクトルとなる．

1.2 1.6 4 0 0.6 0.8 3 0 − −      ₌     ₋  ₋       （42） すなわちアーム軸方向に作用する先端荷重に対しては（図 6 (b)），式（40）の右辺でいかなる駆動力も必要とし ないことを意味する（β ≠0のときの折れたアームでは駆動力は必要となる）．これは軸方向の先端荷重に対し て曲げが発生しないことからも明らかである．以上，双対性のあるJと T J のクリティカルな左右の固有ベクト ルの工学的役割は，図 5 と図 6 において対比するとより明白となる． 6．クリティカルな左特異ベクトルへの拡張 次に力学問題の双対性を維持したまま，正方行列Jに関する議論を長方形行列Bにまで拡張する．ここでも 行列Bの更新に伴い，ランク落ちの原因となる特異値sc(→0)と，s に対応する特異ベクトルを「クリティカc ルな」と形容する．例として図 7 にある 3 本の弾性部材からなる平面トラスを考える．系の自由度は節点 0 の並 進変位( , )U V の 2 個である．節点番号 1, 2, 3 を部材番号 1, 2, 3 として併用する．左右両端の部材 1, 2 の鉛直部材 3 からの開き角ωが徐々に小さくなるつれ，系全体は不安定系に移行するが，当面は式中に微小角ωを残して ( )ω B と表記する．

Fig.7 Statically indeterminate 3-bar truss

2 個の節点変位( , )U V と 3 本の弾性部材の伸び(δ δ δ1, 2, 3)を結びつける幾何学的な適合条件は，次のようである．  1 2 (3 2) 3 ( ) U V δ ω δ δ ×      ₌           B （43） ここで，  (3 2) 1 ( ) 1 0 1 ω ω ω ×     = −_{ }     B （44） であり，有限要素法で節点変位と要素変形（ひずみ）を結びつけるB行列とは本質的には同じである．またB( )ω の特異値分解は，次のようになる．

  _  T 1 T 1 c 0 c (3 1) (3 2) _{(1 2)} (3 3) _{(2 2)} (3 2) 3 0 ( ) 0 2 0 0 ω ω × × _× × _× ×  _{  }    _{  }   =  _{  }      _{ }_ _   v B u u u _v   （45）

ωを残したままのB( )ω のランクはr=2（full column rank）であり，2 個の特異値はs1=sc = 2ωおよびs2= 3

である．開き角ωの変化（ω →0）に伴うB( )ω のランク低下（s1=sc→ ）に対応する左右の特異ベクトルを，0 それぞれu （c 3 1× ）およびv （c 2 1× ）として，次のように定義する．

(

)

T c (0)= 0 0 u B （46） および c 0 (0) 0 0     =       B v （47） ( )ω B のクリティカルな左右の特異ベクトルを抽出するためにそれぞれ，次のような対称行列B B およびT BBT を T LDL 分解する．  2 2 T T (2 2) (2 2) 1 0 1 0 2 0 2 0 0 1 0 1 0 3 0 3 ω ω × ×       = = =            B B _LDL （48） 0 ω → のとき（図 7 (c)）r=1となり，Dの第 1 番目の対角成分がゼロとなるため，B(0)のクリティカルな右特 異ベクトルv （c 2 1× ）はL−Tの第 1 列として抽出できる(7)(9)(10)． c (2 1) 1 0 ×   =     v （49） 同様に_BBT_{も，次のようにLDL 分解する．} T  2 2 2 2 2 2 2 2 T 2 2 T 2 2 (3 3) (3 3) 2 1 1 1 1 0 0 1 0 0 _{1 1} 1 1 1 1 4 1 1 1 1 1 0 0 0 0 1 2 1 1 1 1 1 0 0 1 1 1 0 0 0 1 2 1 ω ω _{ω ω} ω ω ω ω ω ω ω ω ω × ×  −          + _{+ +}  + −        _{ −   } = − + = =_{   } + +         _{ }    _{ } _{ }   ₊    _{ } BB LDL_ （50） T BB についてはω →0のときDの第 2 番目の対角成分がゼロとなるため，B(0)のクリティカルな左特異ベクト ルu はc L−Tの第 2 列として抽出でき(7)(9)(10)， 次のように決まる．

c (3 1) 1 1 0 ×     = −      u （51） 重なる 3 本の部材を分離して描写した図 8 を参照しながら，B(0)のクリティカルな左右の特異ベクトルの工学 的な意味を次に考える．

Fig.8 Critical left and right singular vectors of B Fig.9 Critical left and right singular vectors of B T

0 ω → のとき式（43）から( , )U V を求解するとき，部材の弾性伸び(δ δ δ1, 2, 3)は，任意に設定することはでき ず，クリティカルな左特異ベクトルu と直交する必要がある（図 8 (a)）c ．

(

)

1 2 1 2 3 1 1 0 0 δ δ δ δ δ     − _{ }= → =     （52） すなわち変形の適合条件としてトラス部材の伸びは中央部材 3 に関して左右対称である必要がある（図 8 (a)）． これがB(0)のクリティカルな左特異ベクトルが与える工学的な意味である．ちなみに，u0 =(1,1, 2)− Tとの直交 性からはδ1+δ2 =2δ3が導出され，式（52）と合わせてδ1=δ2 =δ3となる．他方，クリティカルな右特異ベク トル T c=(1, 0) v は，次のような同次解でもある．  c (3 2) 0 (0) 0 0 ×     =       B v （53） したがってv は工学的には，各部材の弾性伸びが生じない剛体変位モードc ( , )U V を意味する（図 8 (b)）． 次に裏の双対問題としての平衡問題に着目する．節点荷重( , )P Q を与えて節点 0 での平衡条件から，未知の部 材力(n1,n2,n3)を求める求解問題は，次のようになる． 1 T 2 (2 3) ₃ ( ) n P n Q n ω ×      _{ =}           B  （54）

この平衡式での T ( )ω B のクリティカルな左右の特異ベクトルは，それぞれv およびc u となる．そしてc ω →0 のとき，クリティカルな左特異ベクトルv と，荷重条件c ( , )P Q との間の直交条件から有意の平衡問題として成立 するためには，P=0が必要とされる．すなわち形態が不安定なトラス系の節点 0 においては，図 9 (a) が示すよ うに，鉛直方向の荷重成分 Q のみが許容できることになる．これが T ( )ω B のクリティカルな左特異ベクトルvc がもつ力学的意義である．ちなみに T ( )ω B のクリティカルな右特異ベクトルu は，次のような平衡方程式の同c 次解でもある． T c (2 3) 0 (0) 0 ×   =     B u  （55） この同次解u は，節点 0 でいかなる外力も作用しないとき，自己平衡状態にあるべき部材内力を与えてくれるc （図 9 (b)）．これが T ( )ω B のクリティカルな右特異ベクトルu がもつ力学的意義である． c このようにB( )ω とBT( )ω との間に力学的な双対性があるとき，図 8 と図 9 の対比からもわかるように，クリ ティカルな左右の特異ベクトルについて，表裏一体化した工学的な意味づけが鮮明となる． 7．応 用（非関連流動則を用いた弾塑性分岐解析） 最後により工学的な応用例として材料の有限要素分岐解析を掲げる． 7・1 弾塑性構成則の概要 数多くの材料パラメータの詳細は，紙面の関係から参考文献(21)~(28)_{に譲ることにして，弾塑性構成則の概要のみ} を示す．対象とした材料は非関連流動則に従う地盤材料であり，材料モデルは，拘束圧依存性やダイレイタンシ ー特性を示す地盤材料を対象として，Drucker–Pragerモデルをベースとして，塑性せん断ひずみによるせん断抵抗 係数の動員，Roweのストレス－ダイレイタンシー則，応力の第三不変量の影響（Lode角の影響）を導入して拡張 されたモデルである．ただし，異方性は考慮していない等方硬化モデルである．弾性部分の構成式には，地盤材 料に適した超弾性モデルを適用した．弾塑性構成式では非関連形の塑性流動則を導入しているため，Kirchhoff応 力のTruesdell速度と変形速度テンソルとの関係として表した速度形構成式の弾塑性接線係数テンソル（4 階テンソ ル）は非対称テンソルとなる(24)~(26)_{．したがって有限要素解析における接線剛性行列も非対称となる}(27)(28)_． 7・2 材料分岐モデルと計算条件 形状比H0 W0 =2.25の矩形平面ひずみ供試体の圧縮解析モデルを図 10 に示す．上下端面は摩擦がない境界と し，上下それぞれの中央節点のみ水平方向変位を拘束し，それ以外の節点では水平方向に自由に変位できる．上 端面が平面を保持する拘束を課して，上端面に圧縮荷重Fを作用させる．このとき平衡路の主経路に沿う変形は 一様変形モードとなる．なお，左右側面に対する側方拘束圧はないものと仮定した．弾塑性分岐解析の手法は， 文献(27)(28)_{と同様であり，詳細は省略する．} 7・3 解析結果（クリティカルな左右の固有ベクトルの検討） 解析の結果，主経路上に計 6 個の特異点が検出された．各特異点において荷重極限点と分岐点の判別を行った 結果，最初の特異点は分岐点，2 番目の特異点は荷重極大点，それ以降の 4 個は分岐点である．このうち荷重極 大点付近に着目し，その拡大図を図 11 に示す．圧縮変位∆Hと圧縮荷重Fをそれぞれ初期寸法H , 0 W で除し0 た公称軸圧縮ひずみと公称軸圧縮応力との関係で示してある．図 11 の荷重極大点（• ）の前後各 1 個の分岐点（ ） をそれぞれ，分岐点（1），（2）と呼ぶことにする．荷重極大点付近で複数の分岐点が出現するため，解析難易度 の高い Hill-Top Branching（頂上分岐）に似た状況にある例題となっている． 分岐点（  ）においてLDU分解された接線剛性行列から抽出されたクリティカルな左右の固有ベクトルを図 12 に示す．図中の太い破線は特異点が検知された時点の変形した供試体形状であり，この変形形状からの変位増

分として左右の固有ベクトルをプロットしている．適切なスケーリングを施しているが，各特異点で左右の固有 ベクトルに対するスケーリング係数は同一とした．右固有ベクトルが分岐座屈モードに対応する．荷重極大点に おけるクリティカルな左右の固有ベクトルが図 13 である．

Fig. 10 Plane-strain compression analysis Fig. 11 Bifurcation points and limit point

図 10 で圧縮荷重Fは供試体の上端面のみに作用していることを考慮すると，図 12 の 2 つの分岐点でのクリテ ィカルな左固有ベクトルと荷重モードベクトルとの間の直交性が確認でき，分岐の瞬間には荷重モードベクトル にクリティカルな左固有ベクトルをあてがった場合に，いかなる外部仕事も生まれないことを意味する．この結 果，可解条件である式（13）が満足されるため分岐解が存在する．図 13 に示した荷重極大点では式（13）のよう な直交性はなく，荷重極大点でちょうど荷重変数の増分量がゼロであることにより，式（12）の右辺におけるb（荷 重ベクトル）はゼロベクトルとなる．このときクリティカルな右固有ベクトル以外に有意な解は存在しないこと は明らかである．なお接線剛性行列の非対称性により各特異点でのクリティカルな左右の固有ベクトルは，程度 の差はあるものの異なるモードとなっている．特に分岐点（2）での差は顕著である．図 12 に示すクリティカル な右固有ベクトルと荷重モードベクトルの間にも直交性がみられるが，この例題での偶然の直交性であり，この ような直交性そのものに意味はない．分岐点と荷重極大点の識別にクリティカルな右固有ベクトルを参照するこ とは，特異点の分類の誤認につながる危険性がある．

(a) Bifurcation point (1) (b) Bifurcation point (2) Fig. 12 Critical left and right eigenvectors at bifurcation points

分岐点（1），（2）から枝分かれする分岐経路も図 11 に示す．分岐経路への切替えを行う際の変位予測子として

各分岐点でのクリティカルな右固有ベクトルを用いた．Clustered 分岐のように複数の分岐点が非常に近接してお り，分岐経路への切替え（Branch-Switching）が困難な条件下であるにも関わらず，各分岐点での分岐モードがク

リティカルな右固有ベクトルとして精確に抽出できていることにより，安定して分岐経路への切替えが可能であ った．本章の例のように，平衡路の追跡過程で荷重極限点と分岐点が混在・近接して発生する問題において効率 的な分岐解析を行うためには特異点の識別と，分岐点である場合には分岐経路への経路切替えのための良質な変 位予測子が不可欠であり，クリティカルな左右の固有ベクトルがそれぞれ重要な役割を果たす．

Fig. 13 Critical left and right eigenvectors at load limit point

8. 結 語 一般に非線形解析において，線形化された方程式の係数行列が特異性を示す工学問題は多い．このような場面 では単に解けないとするのではなく，特異行列の数理構造について考察し，可解にむけての状況を積極的に作り 出すことで，発展的な数理技術の開発ヒントにつながることがある．この際，特にクリティカルな左固有ベクト ルは現実的な工学問題として成立するかどうかを決める重要な役割を演ずる． 工学的応用の例として 7 章でも例示したように，固体や構造の分岐座屈問題(1)(20)_{がある．特異座屈点（または} その近傍）における接線剛性行列は問題により非対称特異行列となる場合もある(11)(12)_{．このとき解くべき剛性方} 程式のひとつは，右辺項に単位荷重増分に対する荷重モードベクトルをもつ．そしてこのような安定問題の解析 においては，飛移り座屈（荷重変数の折り返し点）と分岐座屈（平衡経路の枝分かれ）との数理計算的な識別に 不可欠となるのが，特異剛性行列のクリティカルな左固有ベクトルと荷重ベクトルとの直交条件である．この識 別の際に，クリティカルな右固有ベクトルを用いると，飛移り座屈と分岐座屈とを誤認識する危険性がある． またクリティカルな左右の特異ベクトルの応用例としては，平面応力問題における三角形一定ひずみの CST 要 素（6 自由度）の定式化のなかで，要素内のひずみ 3 成分と節点自由度とを結びつけるB（3 6× ）行列がある． そして_{B （}T 6 3× ）は要素内の応力を節点力に結びつける長方形行列となる．B（3 6× ）のゼロ特異値に対応す る右特異ベクトルは，いかなるひずみをも生まない要素の剛体変位を記述する節点自由度の組み合わせを意味し， ゼロ特異値に対応する左特異ベクトルは，いかなる節点力も生まない自己平衡状態にある応力を意味することに なる． 本研究は計算分岐力学の研究に端を発し，他の多くの工学の数理問題にも共通する本質に着目した基礎研究で ある．なお本文中で提示したすべてのクリティカルな固有ベクトルは固有値解析に依らずに，著者らが既報(7)~(12) の中で提案するように，係数行列の_{LDL 分解または}T LDU分解から抽出している． 文 献 (1) 藤井文夫，大崎 純，池田清宏，構造と材料の分岐力学，計算シリーズ 3 (2005) , p. 192, コロナ社． (2) 川口健一，一般逆行列と構造工学への応用，計算工学シリーズ 1 (2011), p. 214, コロナ社． (3) 藤井文夫，瀧 諭，萩原伸幸，本間俊夫，三井和男，非線形構造モデルの動的応答と安定性，計算工学シリーズ 2 (2003), コロナ社． (4) 京谷孝史，よくわかる連続体力学ノート，非線形 CAE 協会 編 (2008), p. 296 , 森北出版． (5) 清水昭比古，連続体力学の話法 ―流体力学，材料力学の前に― (2012), p. 320, 森北出版．

(6) 野口裕久，久田俊明，“Scaled Corrector を用いた有限要素分岐解析法の開発”，日本機械学会論文集 A 編, Vol. 58, No. 555 (1992), pp. 2191–2198.

(7) 藤井文夫，野口裕久，“固有値解析を必要としない大規模非線形構造系の分岐座屈モードの求め方”，構造工学論文

集（計算力学部門），Vol. 46A (2000), pp. 241–250.

(8) Fujii, F., Ikeda, K., Noguchi, H. and Okazawa, S., “Modified Stiffness Iteration to Pinpoint Multiple Bifurcation Points”,

Computer Methods in Applied Mechanics and Engineering, Vol. 190, Issues 1–9 (2001), pp. 2499–2522.

(9) Fujii, F. and Noguchi, H., “The Buckling Mode Extracted from the LDLT-decomposed Large-order Stiffness Matrix”,

Communications in Numerical Methods in Engineering, Vol. 18, Issue 7 (2002), pp. 459–467.

(10) Noguchi, H. and Fujii, F., “Eigenvector-free Indicator, Pinpointing and Branch-switching for Bifurcation”, Communications in

Numerical Methods in Engineering, Vol. 19, Issue 6 (2003), pp. 445–457.

(11) Fujii, F. and Yamakawa, Y., “Left and Right Eigenvectors Extracted from the LDU-decomposed Jacobian Matrix in Stability Problems”, APCOM2007, Kyoto, CD-ROM Proceedings (2007).

(12) Fujii, F., Yamakawa, Y. and Noguchi, H., “Extracting the Left and Right Critical Eigenvectors from the LDU-decomposed Non-symmetric Jacobian Matrix in Stability Problems”, Computational Mechanics, Springer, Vol. 46, No. 2 (2010), pp. 215–228. (13) 伊理正夫，阪田省二郎，応用数学 3 ＝マトリックス，工学基礎講座 3 (1982), p. 250, 培風館． (14) Strang, G. 著，山口晶哉，井上昭 訳，線形代数とその応用 (1978), p. 440, 産業図書． (15) 木村英紀，線形代数 数理科学の基礎 (2003), p. 248, 東京大学出版会． (16) 柳井晴夫，竹内 啓，射影行列・一般逆行列・特異値分解，UP 応用数学選書 10 (1993), p. 214, 東京大学出版会． (17) 笠原晧司，線型代数と固有値問題 ―スペクトル分解を中心に―，改訂増補版 (2005), p. 332, 現代数学社． (18) 池辺八洲彦，池辺淑子，浅井信吉，宮崎佳典，現代線形代数 ―分解定理を中心に―，初版 (2009), p. 367, 共立出 版． (19) 太田快人，システム制御のための数学（1）―線形代数編―，システム制御工学シリーズ 7 (2000), p. 266, コロナ社． (20) Kobayashi, T., Mihara, Y. and Fujii, F., “Path-tracing Analysis for Post-buckling Process of Elastic Cylindrical Shells under

Axial Compression”, Thin-Walled Structures, Vol. 61 (2012), pp. 180–187.

(21) Yamakawa, Y., Hashiguchi, K. and Ikeda, K., “Implicit Stress-update Algorithm for Isotropic Cam-clay Model Based on the Subloading Surface Concept at Finite Strains”, International Journal of Plasticity, Vol. 26 (2010), pp. 634–658.

(22) 山川優樹，山口洋介，橋口公一，池田清宏，“拡張下負荷面 Cam-clay モデルの有限変形理論に基づく定式化とリタ

ーンマッピングを用いた陰的応力更新法”，応用力学論文集，土木学会，Vol. 13 (2010), pp. 411–422.

(23) Bardet, J.P., “Lode Dependences for Isotropic Pressure-sensitive Elastoplastic Materials”, Transactions of ASME, Vol. 57 (1990), pp. 498–506.

(24) Miehe, C., “On the Representation of Prandtl–Reuss Tensors within the Framework of Multiplicative Elastoplasticity”,

International Journal of Plasticity, Vol. 10, No. 6 (1994), pp. 609–621.

(25) Borja, R.I., “Bifurcation of Elastoplastic Solids to Shear Band Mode at Finite Strain”, Computer Methods in Applied

Mechanics and Engineering, Vol. 191 (2002), pp. 5287–5314.

(26) Hashiguchi, K. and Yamakawa, Y., Introduction to Finite Strain Theory for Continuum Elasto-Plasticity, Wiley Series in Computational Mechanics (2012), John Wiley & Sons.

(27) Ikeda, K., Yamakawa, Y. and Tsutsumi, S., “Simulation and Interpretation of Diffuse Mode Bifurcation of Elastoplastic Solids”,

Journal of the Mechanics and Physics of Solids, Vol. 51, No. 9 (2003), pp. 1649–1673.

(28) Ikeda, K., Yamakawa, Y., Desrues, J. and Murota, K., “Bifurcations to Diversify Geometrical Patterns of Shear Bands on Granular Material”, Physical Review Letters, Vol. 100, No. 19 (2008), Paper No. 198001.