線形システムの安定性を調べる
octaveでisstable( )する話.
追記:octave-forgeのcontrolコンポーネントの中にisstable()関数がありました.そちらを使うことを推奨します.
線形システムの内部安定性
フライトシミュレーションの記事でも扱ったけれど,線形なシステムは一般に
という形で表記できる.特にのときには
となるわけだけど,この時にシステムが安定かどうか(つまり内部安定性)は行列によって決定付けられる.手計算でやるならラウスやフルビッツの方法を使うけれど,計算機があるならそんな手間は不要.そのまま特性方程式の根を求めてやって,根の実部が全部負だったら安定ということになる.
そういった処理を行うとき,Matlabだとisstable( )関数で一発なんだけど,octaveには残念ながらisstable( )関数が無い.poly( )関数で特性多項式を出して,roots( )関数で根を求めて,根の一覧を眺めてナルホドしてもいいんだけど,ちょっと面倒くさい.
octaveでisstable( )する
というわけでisstable( )関数を書いた.行列の特性方程式の根の実部がすべて負であれば安定なので,それを書いてやるだけ.安定ならば1, 不安定ならば0を返す関数になっている.関数スクリプトなので,ファイル名は"isstable.m"にすること.
関数スクリプト
% isstable.m function T_F = isstable(A) %特性多項式を計算 pl = poly(A); %根を計算 r = roots(pl); %安定判別. 根の実部に正のものがあれば不安定. T_F = (length(find(real(r)>0))==0); if (T_F==true) printf('the given system is stable.\n') else printf('the given system is unstable.\n') end end
使い方
行列をisstable( )関数に突っ込む.
octave:1> A = [-1,2;-4,-5]; octave:2> isstable(A); the given system is stable.
安定らしい.確認してみる.
octave:3> roots(poly(A)) ans = -3.0000 + 2.0000i -3.0000 - 2.0000i
実部は確かに負だ.
以上.
Octaveの精義―フリーの高機能数値計算ツールを使いこなす
- 作者: 松田七美男
- 出版社/メーカー: カットシステム
- 発売日: 2010/12
- メディア: 単行本
- クリック: 2回
- この商品を含むブログ (1件) を見る
- 作者: 杉江俊治,藤田政之
- 出版社/メーカー: コロナ社
- 発売日: 1999/01
- メディア: 単行本
- クリック: 13回
- この商品を含むブログを見る