SCILABの基礎事項 | 制御系CAD

[0] はじめに
SCILABは対話型数値計算プログラムで、次からフリーで入手できます。
http://www.scilab.org/products/scilab/download

[1] 行列データの入力
１）サイズ２×２の行列
$A=\left[\begin{array}{cc}1&2\\3&4\end{array}\right]$
を定義するために、次を入力してみます。
–> A=[1 2;3 4]
ここで、行列データは[ ]で囲み、要素の区切りはスペースまたはカンマ、行替えはセミコロン；を用いることに注意してください。次のように最後に；を置くと、出力が抑制されます。
–> A=[1 2;3 4];
すでに入力した行列データへのアクセスは、次のように行います。
–> A
また、特定の要素へのアクセスや修正は、次のように行います。
–> A(2,1)
–> A(2,1)=-5

２）２次の単位行列
$B=\left[\begin{array}{cc}1&0\\0&1\end{array}\right]$
は、次により入力できます。
–> B=eye(2,2)
便利なのは、次のように異なる行数と列数が許されることです。
–> B=eye(2,3)
これは、次のサイズ２×３の行列
$B=\left[\begin{array}{ccc}1&0&0\\0&1&0\end{array}\right]$
の定義を意味します。

３）サイズ３×２の零行列
$C=\left[\begin{array}{cc}0&0\\0&0\\0&0\end{array}\right]$
は、次により定義できます。
–> C=zeros(3,2)

４）３次の対角行列
$D=\left[\begin{array}{ccc}1&0&0\\0&2&0\\0&0&3\end{array}\right]$
は、次により定義できます。
–> D=diag([1 2 3])
すなわち、ベクトルにdiagを作用させているわけですが、行列にdiagを作用させ、対角成分を抽出することもできます。
–> d=diag(D)

５）行列A,B,C,Dを小行列としてもつブロック行列
$E=\left[\begin{array}{cc}A&B\\C&D\end{array}\right]$
を定義するためには、次のように入力します。
–> E=[A B;C D]

６）サイズ２×２の複素行列
$F=\left[\begin{array}{cc}1+i&2+3i\\3+2i&4+4i\end{array}\right]$
も、次により定義できます。
–> F=[1+%i, 2+3*%i;3+2*%i, 4+4*%i]
ここで、%iは虚数単位を意味します。

７）SCILABでは、空行列を定義できます。
–> G=[]

８）すべての要素が1である行列
$H=\left[\begin{array}{cc}1&1\\1&1\end{array}\right]$
は、次により定義できます。
–> H=ones(2,2)

[2] 行列データの管理
１）現時点でどのような行列を定義したかを調べるためには、次のように入力します。
–> who_user
また、たとえば入力した行列データdを消去するためには、次のように入力します。
–> clear d
すべてのデータを消去するときは、次のように入力します。
–> clear

２）行列データEのサイズを調べるには、次のように入力します。
–> [m,n]=size(E)

３）浮動小数点表示を行うためには、次のように入力します。
–> format(’e’)
これを固定小数点表示に戻すためには、次のように入力します。
–> format(’v’)

４）カレントディレクトリのチェックを行うためには、次のコマンドが使えます
–> pwd
–> ls
カレントディレクトリの変更のためには、次のコマンドを与えます。
–> cd(’dname’)

５）メモリ上のデータのセーブとロードのためには、次のコマンドを与えます。
–> save(’fname’)
–> load(’fname’)

[3] コロン記法
１）まず、次を入力してみてください。
–> h=1:3
これは、「form 1 to 3」を意味します。すなわちサイズ1×5の行列h=[1 2 3]が生成されます。コロン「:」はtoの意味で使われ、colon notationと呼ばれています。

２）それでは、次はどうでしょう。
–> h=0:0.5:3
これは、「form 1 to 3 by 0.5」を意味します。すなわち、２つのコロンで挟まれた0.5は増分を意味し、サイズ1×5の行列h=[1 1.5 2 2.5 3]が生成されます。ちなみに、h=1:3はh=1:1:3と同じと考えておきましょう。

３）それでは、次はどうでしょう。
–> h=0:-2
これは、h=1:1:-2の意味になりますが、1から-2まで、増分1では到達できないので、空行列となります。

[4] 小行列の抽出
１）コロン記法の利点は，行列の添え字部分に適用できることにあります。たとえば，行列Eの(2,3)要素には
–> E(2,3)
でアクセスできるのですが，この2と3は，スカラーである必要はなく，ベクトルが使えます。たとえば，次を確かめてください。
–> E([1 3],: )
これは第１行と弟２行からなるサイズ2×5の小行列が取り出されることに注意してください。列を指定する「：」はすべての列を意味しています。

２）次はどうでしょう。
–> E(:,2:4 )
ここで、2:4は[2 3 4]というベクトルでしたよね。第2～4列からなるサイズ5×3の小行列が取り出されることに注意してください。行を指定する「：」はすべての行を意味しています。

３）それでは，次はどうでしょう。
–> E([1 3],2:4 )

[5] 行列の分割
１）いろいろな応用では、ある行列を分割する必要が出てきます。これは小行列の抽出の観点から、たとえば次のようにできます。
–> [m,n]=size(E); i=2; j=2;
–> A=E(1:i,1:j);
–> B=E(1:i,j+1:n);
–> C=E(i+1:m,1:j);
–> D=E(i+1:m,j+1:n);
ここで、i=0のときは1:i=[]、i=mのときはi+1:m=[]などに注意してください。

[6] 行列演算
いま、行列Aが次のように定義されていることを確認してください。
$A=\left[\begin{array}{cc}1&2\\-5&4\end{array}\right]$
また、行列Bを新しく定義します。
–> B=eye(A)

１）行列A,Bの和、差、乗算は次のように計算できます。
–> X=A+B
–> X=A-B
–> X=A*B

２）行列Aが正則のとき、その逆行列は
–> X=inv(A)
で計算されます。しかし、これは逆行列の要素が欲しいとき以外は使わないようにしましょう。
多くの場合、逆行列は次の連立一次方程式
$AX=B$
または
$XA=B$
を解くときに用いますが、SCILABでは、次のように計算します。
–> X=A\B
–> X=B/A
それぞれ左からの割り算、右からの割り算と呼ばれます。
これは数値的に安定な数値計算アルゴリズムを得るためは直交行列だけの演算ですませることが必要だからです。

３）行列Aのべき乗と転置共役は次のように計算でされます。
–> X=A^2
–> X=A’

４）それでは、次の演算結果と比較してみましょう。
–> X=A.*B
–> X=A.\B
–> X=B./A
–> X=A.^2
–> X=A.’
すなわち、演算子の前に「ドット．」をおくと、要素ごとの演算やただの転置になります。

５）次に、組み込み関数としては、次が使えます。
sin,cos,tan
asin,acos,atan
sinh,cosh,tanh
asinh,acosh,atanh
sqrt
exp,log,log10
abs
real,imag

６）行列関数としては次が準備されています。
–> X=expm(A)
–> X=norm(A)
–> X=rank(A)
–> X=det(A)
ぞれぞれ、行列指数関数、ノルム、階数、行列式を計算します。

７）また、固有値問題を解くには次の関数を用います。
–> [V,R]=spec(A)
単に固有値だけを求めるのであれば
–> r=spec(A)
でＯＫです。

８）さらに、特異値問題を解くには次の関数を用います。
–> [U,S,V]=svd(A)
単に特異値だけを求めるのであれば
–> s=svd(A)
でＯＫです。

９）多項式 $x+1$ と $x+2$ の掛け算は、次のように行います。
–> r=convol([1 1],[1 2])

１０）いま、次のベクトルを定義します。
–> x=[3 1 5 4 2]
このとき、次はxの要素を大きい順番に並べ替えます。
–> [y,i]=gsort(x)
yがソーティングされた結果、iが元のどの位置にあったかを示すインデックスです。たとえば
–> x（i)
としてみるとyと同じ結果が出ます。

１１）次は、ベクトルの要素和、最大値、最小値を求める関数です。
–> y=sum(x)
–> y=max(x)
–> y=min(x)

[7] 関数のグラフ
１）次の関数のグラフを描いてみます。
$p(x)=\frac{1}{2\pi}e^{-\frac{x^2}{2}}\quad(-3\le x\le3)$
次のコマンドを与えてください。
–> x=-3:0.1:3;
–> y=1/sqrt(2*%pi)*exp(-x.*x/2)；
–> plot(x,y)
ここで，なぜx*xでなく，x.*xとしているか考えてみてください。

２）描画範囲（-3<=x<=3, 0<=ｙ<=0.4）を指定し，グリッドを描画するには，次のようなコマンドを追加します。
–> mtlb_grid
–> mtlb_axis([-3 3 0 0.5])
–> xtitle(’N(0,1)’)
ここで，「mtlb_」はMATLABのコマンド相当という意味です。

３）次の１次系の周波数応答を片対数グラフを描いてみましょう。
$|\hat{g}(j\omega)|=20\log_{10}\frac{1}{\sqrt{1+\omega^2}}\quad(0.1\le x\le10)$
次のコマンドを与えてください。
–> w=logspace(-1,1);
–> ga=20*log10(sqrt(1+w.^2).^(-1));
–> plot2d(w,ga,logflag=’ln’)
–> mtlb_grid

[8] プログラミング
これまでは対話型入力を行ってきましたが，一連のコマンドをエディタを用いてあらかじめ編集しておいて一括処理する，すなわちscilabでのプログラミングについて学びます。まず，白紙を示すアイコンをクリックして，エディタscipadを起動してください。

１） 1から9までの奇数の和を求める次のプログラムを入力し，ファイル名ex1.sceでセーブします。
//ex1.sce
//—–
s=0;
for i=1:2:9, s=s+i; end
disp(s,’sum=’)
//—–
//eof
scipadのexecuteコマンドをクリックして実行します。答えは25のはずです。

２）マシンイプシロン(machine epsilon)と呼ばれる特別な定数を求めてみましょう。次のプログラムをファイル名ex2.sceでセーブします。
//ex2.sce
//—–
s=1;
while 1+s>1, s=s/2; end
s=s*2;
disp(s,’eps=’),
//—–
//eof
scipadのexecuteコマンドをクリックして実行します。epsは次の%epsと等しくなります。
–-> %eps
次を実行してみてください。
–> 1+%eps/2==1
すなわち，
マシンイプシロン%eps=2^{-k}は，1+2^{-k-1}>1が成り立たなくなる最小のkを用いて定義されます。浮動小数点による数値計算では，1+%eps=1となる正数が存在することに注意してください。

３）次は，if文を使った例です。次のプログラムをファイル名ex3.sceでセーブします。
//ex3.sce
//—–
x=input(’x= ‘);
if imag(x)~=0, disp(’complex number’)
elseif x>0, disp(’positive real number’)
elseif x<0, disp(’negative real number’)
else disp(’zero’)
end
//—–
//eof
scipadのexecuteコマンドをクリックして実行します。

４）関数の例として，行列の階数を計算するものを示しましょう。次のプログラムをファイル名rank2.sceでセーブします。
//rank2.sce
//—–
function r=rank2(A)
s=svd(A); tol=max(size(A))*s(1)*%eps; r=sum(s>tol);
endfunction
//—–
//eof
このとき，次を入力してください。
–> exec(‘rank2.sce’)
–> r=rank2([1 2;3 4])

５）次の微分方程式の解を求めることを考えます。
$\dot{x}(t)=-x(t)+\sin(10t)$
次のプログラムをファイル名ex4.sceでセーブします。
//ex4.sce
//—–
function dx=f(t,x), dx=-x+sin(10*t); endfunction
x0=1; t0=0; t=0:0.01:5;
x=ode(x0,t0,t,f);
plot(t,x), mtlb_grid
//—–
//eof
scipadのexecuteコマンドをクリックして実行します。