Scilabを中心としたMATLABクローン即席入門講座

平成15年12月11日(木)

はじめに

ここではMATLABもどきの数値計算システムのScilabとoctaveの基本的な使い方 を説明します。但し、Scilabの方が機能的に充実している為、Scilabを中心に 実例を示しながら説明します。その為、特に例の環境に言及しない場合は全て Scilab上の例です。但し、octaveやMatlabでも出力形式の違いはあるものの、 使い方と結果には違いはない筈です。尚、ここでの"基本的"の意味は行列の 演算、if文やfor文等の制御文と初歩的な関数の構成に必要な内容とし、 処理言語を用いた大きなシステムの構築、Scilabに付属するSCICOSやMETANET 等のツールの使い方は含めません。あくまでも、ScilabやMATLABとOctaveの 違いが余り無く、一纏めに説明が行える範疇に留めます。尚、実際に入力した 結果はscilabの場合は白のテイブル、Octaveの場合はピンクのテイブルと 色分けを行っています。これは個人的な趣味 (Octave -> Octopus -> 蛸 ->ピンク)の為で、色自体にそれ以上の意味は ありません。

ここでOctaveでファイルを扱う場合には Octaveでfileを利用する話を参照して下さい。又、ScilabのSCICOS に関して、その安直な使い方は SCICOSで遊ぼうを参照して下さい。

次にScilab,octaveを用いる環境はUNIXを前提とします。この即席講座では UNIX以外の環境のScilabでも、外部命令を利用するunix命令を除いて 特に問題は無いと思います。

尚、自分の独習メモ程度の代物の為、"出来る"と記述した個所は確認をして いますが、"出来無い"と記述した個所に関しては、私の理解不足も十分考え られるので、過剰に信頼しないで下さい。この講座に関する批評、問題点や 誤り等の指摘は歓迎します。

  1. Help

    MATLAB,Octave,Scilab共にオンラインヘルプを持っている。共に、help 命令で命令の詳細を調べる事が可能である。Octaveの場合は更に "help -i 事項"で、調べたい事項に関連する事を調べられる様になって いる。例えば、常微分方程式を解きたいがどの様な命令が使えるのか 分らない場合、"help -i ode"と入力すればOctaveのオンライン マニュアルから常微分方程式に関する事項を調べる事が可能である。 尚、octaveではhelpでlessを用いてM-fileの先頭のコメントを表示させ、 help -iでtexinfoのinfo fileを表示させている。その為、octaveを利用 する場合には、これらのGNUのアプリケーションを揃えておく必要がある。 次にScilabの場合はOctave風のオンラインマニュアルではないが、 Scilabウインドウの右端にあるHelpメニューからHelpブラウザ を起動させてキーワード検索を行う事が可能である。この点ではMATLABや OctaveよりもGUIを駆使したものとなっている。

    もし、MATLABの入門書があれば、試しにMATLABの命令をHelpで打ち込ん で見ると面白いだろう。基本的なものでは共通するものが案外多い事や 違った面もより明確になるだろう。

  2. Scilabでの入力

    MATLAB,Octave,Scilabの四則演算は通常の数学記号で表現される。 入力行末に記号";"をつけておくと計算結果のエコーバックを行わない。 これは大きな行列データ処理や反復回数が多い処理で必須となるだろう。

    -->1+1
     ans  =
     
        2.  
     
    -->1+1;
     
    -->1+1:
    
    -->
    		

    代入は=で行う。一部の数式処理の様に:=で代入を行わないので混同を しない事。

    -->a = 1.2  
     a  =
     
        1.2  
     
    -->a        
     a  =
     
        1.2  
     
    -->b = %pi*3
     b  =
     
        9.424778  
     
    -->c=a-b
     c  =
     
      - 8.224778  
    
    -->
    		
  3. 数学的定数

    円周率π、自然底e及び純虚数i等の数学的定数はMATLAB,Octave,Scilabに 用意されているが、その値は立ち上げ時に初期化ファイルを用いて設定 されるものである。尚、MATLABやOctaveではそれらの変数名(実際は組込 関数)はpi,e,i等であるが、Scilabでは特にこの様な定数には頭に%を付け、 %pi,%e,%i等の様に通常の変数と予め区別している。その為、FOR文でi=1等 と間違って代入して、その後で純虚数を再定義する事が生じ難くなって いる。OctaveとMATLABでは%が註釈行の先頭とみなされる為、この方法は 使えない。因に、ScilabではJavaと同様に先頭が//で始まる行が註釈行と 見なされる。以下にScilabでの計算例を示す。

    -->%pi 
     %pi  =
     
        3.1415927  
     
    -->%i
     %i  =
     
        i    
     
    -->%e
     %e  =
     
        2.7182818  
     
    -->%i*%i
     ans  =
     
      - 1.  
    
    -->a = 1 + %i;
     
    -->b= 1-%i;a*b
     ans  =
     
        2.  
     
    		

    次にOctaveで同様の計算を行った場合を示す;

    octave:1> pi
    pi = 3.1416
    octave:2> i
    i = 0 + 1i
    octave:3> e
    e = 2.7183
    octave:4> i*i
    ans = -1
    octave:5> a=1+i;
    octave:6> b=1-i;a*b
    ans = 2
    octave:7> 
    		

    前述の様にOctaveではi,e,piは予約語として特別に保護されていない事に 注意が必要である。例えば、Octaveで変数iに適当な値を入れる例を以下 に示す;

    octave:1> who
    octave:2> type i
    i is a builtin function
    octave:3> i=2
    i = 2
    octave:4> i*i
    ans = 4
    octave:5> type i
    i is a user-defined variable
    2
    octave:6> type pi
    pi is a builtin function
    octave:7> pi=10
    pi = 10
    octave:8> i=sqrt(-1)
    i = 0 + 1i
    octave:9> pi=acos(-1)
    pi = 3.1416
    octave:10> who
    
    *** local user variables:
    
    i   pi
    
    octave:11> type pi
    pi is a user-defined variable
    3.1416
    octave:12> for i=1:10
    > i*2;
    > end;
    octave:13> i
    i = 10
    
    octave:14> type i
    i is a user-defined variable
    10
    		

    この様にOctaveでは無警告で変更されてしまうので注意が必要である。 尚、一番最初にwhoを実行した時点では何も表示されなかったが、iとpi に値を入れると利用者定義の変数になるので、whoを実行するとiとpiが 表示されている。更に、typeで調べると利用者定義変数であると表示さ れている。次のforの例では、iが最終的に10で置換えられている事を示し ている。この様にiはindexしてfor文等で迂濶に利用する可能性が高い為、 変数の末尾に1,2,3.. とか数字を付けたり、iiの様に変数は全て二文字 以上にすると云った工夫が必要だろう。

    便利な定数に%epsがある。この定数は微少値(scilabの場合は %eps=4.441E-16程度、ocatveの場合はeps=2.2204e-16とオーダーは同じ。 計算機環境によって異なるだろう)を表現しており、零による割算を避け る為や、while文等のLoopで%epsの何倍かに誤差が収まった場合にLoopを 抜ける様に設定する等と何かと便利な定数である。

  4. 四則演算

    MATLAB、octave、Scilab等での四則演算は通常のプログラム言語と変わ らない形で入力を行う。

    記号 演算
    + 1 + 2
    - 1 - 2
    * 2 * 4
    / 2 / 3
    ^ 羃乗 3^2

    次に、注意する点としてはMATLAB系の言語はCと比べて処理速度が遅い事 である。これは行列の演算のみに関しては目立つものでは無いが、 FOR文によるループを多用した場合には顕著になる。ここで、四則演算 に関しても対処方法を誤れば、計算時間を無駄に消費し兼ねない。 因に、上記の演算の順が丁度、処理速度の速さの順番となり、羃乗は他の 演算と比べて極端に処理が遅い。例えば、a^2やa^0.5はa.*aとsqrt(a)の 方が圧倒的に高速である。この例をOctaveで示す。この例ではrandで 1000行1列の行列を生成し、a.*a,a.^2で成分の二乗を計算している。 尚、返されている数値が処理時間である。尚、この計算は Pentium!!! 600MHzのLinux PCで実行した結果である。但し、この数値 は絶対的なものでは無く、参考値として挙げておく。

    octave:1> a=rand(1000,1);
    octave:2>  t1=time;a.*a;t2=time;t2-t1
    ans = 0.00055897
    octave:3>  t1=time;a.^2;t2=time;t2-t1
    ans = 0.0010910
    octave:4>  t1=time;exp(2*log(a));t2=time;t2-t1
    ans = 0.0024740
    octave:5>  t1=time;a.*a.*a.*a.*a.*a.*a.*a.*a.*a.*a.*a;t2=time;t2-t1
    ans = 0.0011801
    octave:6>  t1=time;a.^12;t2=time;t2-t1
    ans = 0.0010090
    octave:7>  t1=time;exp(12*log(a));t2=time;t2-t1
    ans = 0.0022840
    		

    この様に、二乗の計算でさえも、4倍近くの違いで通常の積の方が 羃よりも速い。又、指数関数と対数関数の組み合わせのものと比べる と通常の積は5倍近い。而し、羃の指数が増えても羃や対数関数を用い たもので差は、処理時間に大きな差が生じないが、積に関しては、 回数が単純に増加する為に、計算時間に大きな違いが出て、羃よりも 遅くなる。

    これはScilabでも同様の傾向を示す。この例でもrandで1000行1列の ベクトルに対して積と羃乗による差の計算を示している。ここでは for文によるloopを利用している。この理由はloop無しの場合、timer() で計測すると0.を返した為である。尚、ここでの計測は Pentium!!! 600MHzのLinux環境での結果である。

    -->a=rand(1000,1);     
     
    -->timer();for k=[1:1000]; b=a.*a;end;timer()     
     ans  =
     
        0.06  
     
    -->timer();for k=[1:1000]; b=a.^2;end;timer()
     ans  =
     
        0.51  
     
    -->timer();for k=[1:1000]; b=exp(2*log(a));end;timer() 
     ans  =
     
        1.26  
     
    -->timer();for k=[1:1000]; b=a.*a.*a.*a.*a.*a.*a.*a.*a.*a.*a.*a;end;t
     ans  = 
     
        0.36  
     
    -->timer();for k=[1:1000]; b=a.^12;end;timer()
     ans  =
     
        0.61  
     
    -->timer();for k=[1:1000]; b=exp(12*log(a));end;timer()
     ans  =
     
        1.3  
    		

    この様に羃乗の処理は遅い傾向を持つ。但し、羃の指数が大きくなった 場合には、安易に積を計算すると積の個数に比例して計算時間が増加す ると考えられるが、羃に関しては極端な差が出る訳では無い。その一方 で、二乗で計算を纏める等の工夫を行う事で全体の計算量の増加を押え る事も可能である。この例をOctaveとScilabで示す。

    先ず、Octaveの例を以下に示す。この例では与えられたリストの 各成分の9乗を計算するものである。

    octave:118> t1=time;b=a.*a.*a;b=b.*b.*b;t2=time;t2-t1
    ans = 0.0057989
    octave:119> t1=time;b=a.^9;t2=time;t2-t1
    ans = 0.0082800
    octave:120> t1=time;b=a.*a.*a.*a.*a.*a.*a.*a.*a;t2=time;t2-t1
    ans = 0.0091070
    octave:121> t1=time;exp(9*log(a));t2=time;t2-t1
    ans = 0.014720
    octave:122> 
    		

    この結果では、羃と安易な積で羃の方が僅かに速くなっている。 これに対し、三乗で纏めた積は羃よりも1.5倍程度高速である。 この様にOctaveでは工夫次第で積による処理の方が速い事が分る。 次に、Scilabで同じ計算を行った例を示す。

    
    -->timer();for k=[1:1000]; b=a.*a.*a;b=b.*b.*b;end;timer()
     ans  =
     
        2.76  
     
    -->timer();for k=[1:1000]; b=a.^9;end;timer()             
     ans  =
     
        2.41  
     
    -->timer();for k=[1:1000]; b=a.*a.*a.*a.*a.*a.*a.*a.*a;end;timer()
     ans  =
     
        3.75  
     
    -->timer();for k=[1:1000]; b=exp(9*log(a));end;timer()            
     ans  =
     
        4.96  
     
    		

    この計算結果ではOctaveと異った結果が得られている。この結果からは Scilabでは少し事情が異なり、羃乗による演算の方がやや速い事になる。 しかし、両方共、logを利用した場合は他と比べ速度的に劣る。

    但し、Pentium3 600MHzのLinux環境ではOctaveど同様の以下の結果を 得ている。

    
     -->a=rand(1000,1);
     
    -->timer();for k=[1:1000]; b=a.*a.*a;b=b.*b.*b;end;timer()
     ans  =
     
        0.16  
     
    -->timer();for k=[1:1000]; b=a.^9;end;timer()       
     ans  =
     
        0.41  
     
    -->timer();for k=[1:1000]; b=a.*a.*a.*a.*a.*a.*a.*a.*a;end;timer()
     ans  =
     
        0.2  
     
    -->timer();for k=[1:1000]; b=exp(9*log(a));end;timer()           
     ans  =
     
        0.85  
     
    -->
    
    		

    この例でも、積を工夫したものが一番速い事が分る。これらの例から、 指数が大きくない羃乗とlog関数は可能な限り避けた方が賢明である。 ここで羃乗に関しては指数の大きさによる面もある為、注意が必要だが、 全体的な計算量を削減したものの方が良好な結果を得易い。この様に 全体の計算量を考慮せずに、安易な処理を行った場合には、逆に悪化 する可能性が高くなる。尚、それ以上の処理速度がどうしても必要ならば 別の手段を考える必要がある。例えば、行列の積演算に関してATLASを 用いて高速化する方法、更にはC等で計算の核となる部分を記述し、 リンクして利用する事を検討した方が良い。因に、Scilab、Octave、 MATLAB共に外部のCやFortranプログラムを組み込む事が可能である。 但し、この組み込み方法はScilab、Octave、MATLABで全く異る。 特に、MS-Windows環境上のMATLABであれば、そのコンパイラにも注意 する必要がある。尚、この組み込み方法に関してはここでは述べない。

    又、処理速度に関しては不必要なfor文を用いれば、容易に低下させる事 が可能である。その為、行列の処理等ではfor文を用いた成分毎の処理を 避け、行列同士の演算か、行列をベクトルの列又は行と見なして処理を 行う方が処理速度で優位にある。勿論、非常に大きな行列の直接演算は それ相応の処理時間を必要とし、疎行列を扱う場合は特に、直接演算 よりも、それに適した演算処理を行う方が良い。但し、for文の利用を 出来るだけ避けて、行列やベクトル単位の演算を主体で行う方が処理 速度の面で優位にある。これはMATLAB言語等では、その言語仕様自体が、 先ず、配列から指定された成分を取り出してコピーし、そのコピーを 用いて計算した結果を再び配列に戻す事を行っている。これに対し、 MATLABやScilabでは行列やベクトルの計算でBLAS等の標準的なライブラリ を直接利用する為、配列の取り出しやコピーの様な余計な手間が不要と なる為である。

  5. 行列

    MATLABクローンの言語で扱うデータは基本的に行列である。言語仕様 自体が効率良く行列を扱える様に工夫されている。CやFORTRANの配列は ベクトル又は行列として表現される。ここで注意する事の一つに指数の 開始番号である。Cでは配列番号は0から始まるが、MATLABクローンでは FORTRAN同様に番号は1から開始する。この点は数学的な表現と違いが 無いので混乱する事は無いだろう。又、ベクトルや行列の入力はScilab、 MATLAB、Octaveでも同じ形式で行う。行列の入力例を以下にScilabを 用いて示す。

    -->[1,2,3;4,5,6]
     ans  =
     
    !   1.    2.    3. !
    !   4.    5.    6. !
    
    -->[1;2;3]
     ans  =
     
    !   1. !
    !   2. !
    !   3. !
    
    -->a=[1 2 3]
     a  =
     
    !   1.    2.    3. !
    
    		

    これらの例に示す様にMATLAB,Octave,Scilabでは行列やベクトルに関し てキャラクターを用いた数学的/視覚的な表示を行う事も特徴の一つ である。又、元の区切はspaceか","を用い、行の区切には";"を用いる。 但し、";"に似ている":"は別の意味を持つ。この":"の便利な使い方に 関しては後に述べる。

    次にScilabやMATLABでは、入力した行列aのi行j列の成分をa(i,j)で表現 する(a[i,j]では無い事に注意)。この際、i,jは数学の行列表記と同様に 1から始まり、Cの様にから開始するのでは無い。ベクトルや行列の元の 設定は上記の例の様に一気に設定する方法と、a(1,2)=1の様に成分を 一々指定して設定する方法もある。この場合、aが未定義で、最初に a(1,2)=1とすると、aは1行2列の行列として生成され、1行2列の成分 a(1,2)のみが1で他が零となっている。更に、a(3,2)=10と続けて入力する と、aは3行2列の行列となり、値が代入された個所(a(1,2)とa(3,2))を除い て、値は零となっている。この方法以外に行列をベクトルとみなして 要素を指定する事も可能である。具体的には、aがm行n列の行列であれば、 a(i,j)はa((j-1)*m+i)として解釈される。但し、この場合はaが何等の 方法で予め定義された場合に限る。aが未定義の場合、この方式でaの 各成分に値を設定すると、aは単純にそのままベクトルとなる。 このベクトルが縦か横かは利用しているMATLABクローン環境に依存する。 ここで行列の大きさはsize命令で行列としての大きさが返される。 尚、MATLABやOctaveとScilabでは行列に対するlength命令がやや異なる。 length命令はベクトルに対してはその長さを返すものである。しかし、 行列(但し、列ベクトル、行ベクトルの何れでも無いとする)に対しては 意味が異なり、Octaveの場合は行数を返し、Scilabの場合は行列の 要素数(ベクトルとみなした場合の長さ)を返す。以下の例で、lengthが 行列の総数を返している事に注目されたい。

    -->a=[1 2 3;4 5 6]
     a  =
     
    !   1.    2.    3. !
    !   4.    5.    6. !
     
    -->a(1,2)
     ans  =
     
        2.  
     
    -->a(2,1)
     ans  =
     
        4.  
     
    -->a(4)
     ans  =
     
        5.  
     
    -->a(3)
     ans  =
     
        2.  
     
    -->size(a)
     ans  =
     
    !   2.    3. !
     
    -->length(a)
     ans  =
     
        6.  
    		

    行列の一部を取り出す事はScilabやMATLABでは非常に容易に行える。 先ず、ベクトルに対してi番目の成分からj(>=i)番目の成分を持つ ベクトルはa(i:j)で得られる。同時に、行列Aのi行目で構成される 行ベクトルとj列目で構成される列ベクトルはそれぞれA(i,:)と A(:,j)で得られる。尚、最初に注意した様にCの配列と異なり、 番号は通常の行列と同様に1から始まる。次の例では定義した行列 から列ベクトルや行ベクトルを取り出す事を行っている。

    -->a=[1:10]        
     a  =
     
    !   1.   2.    3.    4.    5.    6.    7.    8.    9.    10. !
     
    -->a(4:6)
     ans  =
     
    !   4.    5.    6. !
    
    -->A=[1:4;5:8;8:11]
     A  =
     
    !   1.    2.    3.     4.  !
    !   5.    6.    7.     8.  !
    !   8.    9.    10.    11. !
    
    -->b=A(2,:)        
     b  =
     
    !   5.    6.    7.    8. !
     
    -->c=A(:,2)
     c  =
     
    !   2. !
    !   6. !
    !   9. !
     
    -->
    		

    ところで、記号":"を用いて行列同士の代入を簡略化する事が可能である。 先ず、先程示した様に行列Aが与えられた場合、A(:,1)が行列Aの一列目の 元で構成される列ベクトルを意味し、A(2,:)は行列Aの二行目の成分で 構成される行ベクトルを意味する。次の例では行列Aの第一列A(:,1)と 第二行A(2,:)を取り出している。

    -->A=rand(4,4)
     A  =
     
    !   0.2113249    0.6653811    0.8782165    0.7263507 !
    !   0.7560439    0.6283918    0.0683740    0.1985144 !
    !   0.0002211    0.8497452    0.5608486    0.5442573 !
    !   0.3303271    0.6857310    0.6623569    0.2320748 !
     
    -->A(:,1)
     ans  =
     
    !   0.2113249 !
    !   0.7560439 !
    !   0.0002211 !
    !   0.3303271 !
     
    -->A(2,:)
     ans  =
     
    !   0.7560439    0.6283918    0.0683740    0.1985144 !
     
    -->
    		

    この記号":"を利用して、行列を列或いは行ベクトルの集まりとみなして、 これらの成分を持つ行列を以下の例の方法で容易に生成出来る。

    -->B=zeros(4,4)    
     B  =
     
    !   0.    0.    0.    0. !
    !   0.    0.    0.    0. !
    !   0.    0.    0.    0. !
    !   0.    0.    0.    0. !
     
    -->B(:,2)=A(:,4)
     B  =
     
    !   0.    0.7263507    0.    0. !
    !   0.    0.1985144    0.    0. !
    !   0.    0.5442573    0.    0. !
    !   0.    0.2320748    0.    0. !
     
    -->
    	      

    但し、この代入で注意する事はa(:,1)で代入を行う場合、代入される 側のサイズが":"で省略されている場合、両方とも":"で示す大きさは 同じものでなければならない。

    -->C=zeros(2,2)
     C  =
     
    !   0.    0. !
    !   0.    0. !
     
    -->C(:,2)=A(:,2)
                  !--error    15 
    submatrix incorrectly defined
     
    -->a(:,1)=A(:,4)
                  !--error    15 
    submatrix incorrectly defined
    
    -->a=A(:,4)
     a  =
     
    !   0.7263507 !
    !   0.1985144 !
    !   0.5442573 !
    !   0.2320748 !
    		

    この上の例に示す様に":"で示した行列が右辺と左辺の大きさが違う場合 にエラーになる。又、左辺が未定義の場合に":"を用いて成分を指定する 場合もエラーになる。この場合には予め変数を定義しておく必要がある。 但し、右辺と左辺の大きさが一致すれば問題が無い。

    以下に例を示す;

    -->D=zeros(10,2)
     D  =
     
    !   0.    0. !
    !   0.    0. !
    !   0.    0. !
    !   0.    0. !
    !   0.    0. !
    !   0.    0. !
    !   0.    0. !
    !   0.    0. !
    !   0.    0. !
    !   0.    0. !
     
    -->A 
     A  =
     
    !   0.2113249    0.6653811    0.8782165    0.7263507 !
    !   0.7560439    0.6283918    0.0683740    0.1985144 !
    !   0.0002211    0.8497452    0.5608486    0.5442573 !
    !   0.3303271    0.6857310    0.6623569    0.2320748 !
     
    -->D([7:10],1)=A(:,1)
     D  =
     
    !   0.           0. !
    !   0.           0. !
    !   0.           0. !
    !   0.           0. !
    !   0.           0. !
    !   0.           0. !
    !   0.2113249    0. !
    !   0.7560439    0. !
    !   0.0002211    0. !
    !   0.3303271    0. !
     
    -->
    		

    この様に両者の大きさが異なる場合に両辺に":"を用いて代入すれば、 エラーになるが、両辺のサイズが同じになる様にindexを調整して いれば代入を行う際に問題が無い。これはMATLABで生じる事だが、 代入する変数の大きさを予め指定せず行う場合、右辺と左辺のサイズが 異る事である。具体的には、bをm行、n列の行列として変数aに直接代入 を行う時、a(:)=b(:,2)の様に代入するとaが右辺と同じ列ベクトルと ならずに行ベクトルとなってしまう。尚、これは":"では無く、より 一般のindexを指定しても、行列としての成分を明記していなければ 生じてしまう。例えば、a([1:10])=b([20:-1:11],2)の様に代入した 場合も左辺が行ベクトルになってしまい、両辺のベクトルの形態が 異なる。但し、a=zeros(m,n)の様にaのサイズを定めた場合や、 a=b(:,2)とa(:,1)=b(:,2)の様に左辺の成分を指定せずに直接代入を 行う場合と左辺の成分も指定した場合は両辺共に列ベクトルとなる。 但し、Scilabでは、いきなりa(:)=b(:,2)の様な代入は出来ないが、 a([1:4])=b(:,2)の様な代入は両辺のサイズが同じであれば可能で、 MATLABの様に勝手にベクトルのタイプが変更される事は無い。

    -->aa=rand(4,3)
     aa  =
     
    !   0.2113249    0.6653811    0.8782165 !
    !   0.7560439    0.6283918    0.0683740 !
    !   0.0002211    0.8497452    0.5608486 !
    !   0.3303271    0.6857310    0.6623569 !
     
    -->c([1:4])=aa(:,2)
     c  =
     
    !   0.6653811 !
    !   0.6283918 !
    !   0.8497452 !
    !   0.6857310 !
     
    -->dd(:)=aa(:,3)
                  !--error    15 
    submatrix incorrectly defined
     
     
    -->b=aa(:,1)
     b  =
     
    !   0.2113249 !
    !   0.7560439 !
    !   0.0002211 !
    !   0.3303271 !
    		

    MATLABで上述の現象が生じる理由には、どうもMATLABでは行ベクトルを 基本的データ列と勝手に考えている節があり、その一因に行ベクトルを 表示するとOctaveの様にlessを用いてスクロールを止める機能が無い事 に関係している様に思える。尚、octaveでは代入の際に行ベクトルを 勝手に列ベクトルに変換して代入する事は行わない。その為、この代入が octaveで動作するプログラムがMATLABで正しく動作しなくなる原因の 一つになる可能性が高い。これを避ける為に、aの大きさを予めzerosや onesで初期化、a=b(:,2)の様な直接代入、a(:,3)=b(:,2)の様に要素を 明示して代入する必要がある。

    この様に記号":"は非常に便利な記号である反面、代入で利用する場合 には両辺のサイズが一致している事や、MATLABでの利用を考えている 場合は、尚更、勝手に列ベクトルに変換される事が無い様に注意して 利用する必要がある。その上、配列全体を表す事も出来れば、 定数i,j(>=i)に対し、i:jと記述する事で、先頭がi、末端をjとする 1刻みの数列が得られる。この刻み幅の指定も可能で、例えば、刻み幅 をdとする場合にはi:d:jとすれば良い。但し、j:iとすると、仮定から iの方がjよりも小さい為に、増分1でjからiに到達出来ずに空行が 返される。この場合は刻み幅を負の数にする。更に、1:2:6の様に 大小関係に問題は無いが刻み幅dでiからjに到達出来ない場合、 jを越えない最大の数で数列が止まる。これらの例を以下に示す。

    -->1:5  
     ans  =
     
    !   1.    2.    3.    4.    5. !
     
    -->1:2:5
     ans  =
     
    !   1.    3.    5. !
    
    -->1.1:0.1:1.5
     ans  =
     
    !   1.1    1.2    1.3    1.4    1.5 !
     
    -->5:1   
     ans  =
     
         []
     
    -->5:-2:1
     ans  =
     
    !   5.    3.    1. !
     
    -->1:2:6 
     ans  =
     
    !   1.    3.    5. !
    
    -->1:0.3:2
     ans  =
     
    !   1.    1.3    1.6    1.9 !
    
    		

    行列の演算も数値の四則演算と同様に行える。但し、"*"と".*"、 及び"/"と"./"は別物である。"."が付いたものは行列同士の項単位の 演算を表現しており、この場合は可換である。行列の積はあくまでも "*"である。但し、行列の割り算は非可換の為、二の記号"\","/"が 用意されているので注意する。これらの記号"/","\"は次で説明する。

    -->a = [ 1,2;3,1]
     a  =
     
    !   1.    2. !
    !   3.    1. !
     
    -->b = [1;0]
     b  =
     
    !   1. !
    !   0. !
     
    -->a*b
     ans  =
     
    !   1. !
    !   3. !
     
    -->c = [ 2,3;1,2];
     
    -->d = a/c
     d  =
     
    !   0.    1. !
    !   5.  - 7. !
     
    -->d*c-a
     ans  =
     
    !   0.    0. !
    !   0.    0. !
    
    		

    次に各元に対する積、商及び羃乗を行う".*"、"./"と".^"の実例を示す;

    -->[1 2 3 ].*[3,2,1]
     ans  =
     
    !   3.    4.    3. !
     
    -->[1 2 3 ]./[3,2,1]
     ans  =
     
    !   0.3333333    1.    3. !
     
    -->[1,2,3;3 2 1]./[1,2,3;3,2,1]
     ans  =
     
    !   1.    1.    1. !
    !   1.    1.    1. !
     
    -->[1,2,3;3 2 1].*[1,2,3;3,2,1]
     ans  =
     
    !   1.    4.    9. !
    !   9.    4.    1. !
    
    -->[1,2,3;3 2 1].^[4,5,6;3,2,1]
     ans  =
     
    !   1.     32.    729. !
    !   27.    4.     1.   !
     
    -->[1,2,3;3 2 1].^2            
     ans  =
     
    !   1.    4.    9. !
    !   9.    4.    1. !
     
    -->
    
    		

    上の二つは行ベクトル同士、下の二つは行列同士の項単位の演算を行って いるが、項単位の演算を行う際に、演算を行うベクトル/行列が共に同じ サイズである事に注意する。これは同様に項単位の演算となる"+"と同様 である。但し、大きな行列から一部を同じサイズ程取り出してこの演算を 行うのは一向に構わない。以下にその例を示す;

    -->a=[1:10];
     
    -->b=[1:10;10:-1:1];
     
    -->a(4:8).*b(1,1:5)
     ans  =
     
    !   4.    10.    18.    28.    40. !
    
    -->a(4:8)./b(1,1:5)
     ans  =
     
    !   4.    2.5    2.    1.75    1.6 !
     
    -->a(4:8).^b(1,1:5)
     ans  =
     
    !   4.    25.    216.    2401.    32768. !
    
    -->a(4:8).^2       
     ans  =
     
    !   16.    25.    36.    49.    64. !
     
    -->1./b(1,1:5)     
     ans  =
     
    !   0.0181818 !
    !   0.0363636 !
    !   0.0545455 !
    !   0.0727273 !
    !   0.0909091 !
     
    		

    この例では行列aから4列目から8列目の5成分の行ベクトル[4,5,6,7,8]と 行列bから1行目の1列目から5列目の成分で構成される行ベクトル [1,2,3,4,5]を取り出し、それらの要素毎の積を計算している。 この様に項単位の演算".*"./",".^"はどちらか一方が定数の場合か、 両方とも同じサイズの行列(部分行列も含む)に対してのみ行える。

    次に、行列Aの逆行列はinv命令でinv(A)とするか、単純にA^(-1)で計算 出来る。逆行列の性質上、正則行列のみ計算可能である。又、行列A、B に対し、A^(-1)*bとA*B^(-1)の計算は各々A\B,A/Bでも処理可能である。 但し、MATLABもどきの処理言語が幾ら行列演算向けだからと言っても、 行列の大きさに無関係に逆行列を安易にA^(-1)で計算する事は薦めない。 行列が大きくなれば、それなりの工夫も必要となる。勿論、この事で 考え込む事になるのは、最近のPentium !!!程度の計算機であれば、 1000x1000以上の行列を扱う場合だろう。又、100x100程度でも、処理の 高速化を計りたい場合も同様である。勿論、この辺の事情は計算機環境 に大きく依存する。その為、rand命令で適当なサイズの行列を生成して 予め計測しておくと良いだろう。

  6. ベクトルに対するパターンマッチング

    ScilabやMATLABではベクトルに対してパターンマッチング処理が可能 である。このパターンマッチングを適用する事によって、for文等の loop文を使わずに、見通しの良いプログラムを行う事が可能になる。 この機能は数値データを扱う場合には非常に強力であり、MATLAB クローンの威力が発揮される機能でもある。

    この機能を利用すれば、与えられたベクトルから適合するものがあるか どうか検証する事も可能である。以下の例では与えられたベクトルから 2に等しいものがあるか検証し、その場所をfind命令を用いて探す処理 を行っている。尚、Scilabでは以下に示す様に、x==2の結果として表わ れるBooleanはTとFの値を持っている。Booleanを数値に変換するとTが1、 Fに0が各々対応する。尚、OctaveとMatlabではTが1、Fが0と直接数値が 返される。その為、後で紹介するパターンマッチングで発生した行列 データを用いた演算処理では、Scilabにデータ型の変更が内部的に行なわ れている事に注意されたい。

    -->x=[1:5,5:-1:1]     
     x  =
     
    !   1.    2.    3.    4.    5.    5.    4.    3.    2.    1. !
     
    -->x==2
     ans  =
     
    ! F T F F F F F F T F !
     
    -->y=find(x==2)
     y  =
     
    !   2.    9. !
     
    -->x(y)
     ans  =
     
    !   2.    2. !
     
    -->
    		

    上記の例では2と等しいものを探したが、Cと比べても非常に簡単な処理で 2と等しい元の位置を見付け出している。この検出は等号だけでは無く、 不等号に対しても用いる事が可能である。

    -->x=[1:5,5:-1:1]     
     x  =
     
    !   1.    2.    3.    4.    5.    5.    4.    3.    2.    1. !
     
    -->x>3                
     ans  =
     
    ! F F F T T T T F F F !
     
    -->y=find(x>3)
     y  =
     
    !   4.    5.    6.    7. !
     
    -->x(y)       
     ans  =
     
    !   4.    5.    5.    4. !
    
    -->y=x(x>3)
     y  =
     
    !   4.    5.    5.    4. !
     
    -->x.*y
     ans  =
     
    !   0.    0.    0.    4.    5.    5.    4.    0.    0.    0. !
    
    		

    このfind命令は与えられた行列で0と異なる成分の位置を返す命令である。 但し、Octaveで行列は(i,j)で指定しなければならない為、パターン マッチングを行うものがベクトルでなければx(find(x>3))の様な処理は エラーになる。尚、find(x>3)で返された列ベクトルはx>3で生成 された行列を列ベクトルから構成されたものとみなしている。具体的には mxn行列の(i,j)成分はm*n成分の列ベクトルのm*(j-1)+i番目の成分に対応 する。

    octave:5> aa=rand(5);
    octave:6> bb=aa>0.5 
    bb =
    
      1  1  0  0  1
      1  1  0  1  0
      0  1  1  0  1
      1  0  0  1  0
      0  1  1  1  0
    
    octave:7> find(bb)
    ans =
    
       1
       2
       4
       6
       7
       8
      10
      13
      15
      17
      19
      20
      21
      23
    
    octave:8>aa(find(bb))
    error: single index only valid for row or column vector
    error: evaluating index expression near line 8, column 1
    octave:8> 
    
    		

    而し、Scilabの場合は行列に対してもaa(find(bb))の様な事が可能だが、 答は列ベクトルとなっている。この計算結果はOctaveと同じ配列の対応 となっている。

    -->aa=rand(5,5);
     
    -->bb=aa>0.5
     bb  =
     
    ! F T T F F !
    ! T T T F T !
    ! F T T F F !
    ! F T F T F !
    ! T F T T F !
     
    -->aa(bb)
     ans  =
     
    !   0.7560439 !
    !   0.6653811 !
    !   0.6283918 !
    !   0.8497452 !
    !   0.6857310 !
    !   0.8782165 !
    !   0.5608486 !
    !   0.6623569 !
    !   0.7263507 !
    !   0.5442573 !
    !   0.8833888 !
    !   0.6525135 !
    !   0.9329616 !
     
    -->aa
     aa  =
     
    !   0.2113249    0.6283918    0.5608486    0.2320748    0.3076091 !
    !   0.7560439    0.8497452    0.6623569    0.2312237    0.9329616 !
    !   0.0002211    0.6857310    0.7263507    0.2164633    0.2146008 !
    !   0.3303271    0.8782165    0.1985144    0.8833888    0.312642  !
    !   0.6653811    0.0683740    0.5442573    0.6525135    0.3616361 !
    		

    MATLABクローンでは、y=x(x>3)の様にfindを利用せずに処理する事も 可能である。この様にこのパターンマッチングを適用して処理の簡略化が 行える。例えば、上記の与えられたベクトルから3よりも大きな数値に対し てのみ2倍すると言った計算が次の様に行える。

    -->x=[1:5,5:-1:1]     
     x  =
     
    !   1.    2.    3.    4.    5.    5.    4.    3.    2.    1. !
     
    -->y=find(x>3)
     y  =
     
    !   4.    5.    6.    7. !
    
    -->for i=x(y)
    -->2*i
    -->end
     ans  =
     
        8.  
     ans  =
     
        10.  
     ans  =
     
        10.  
     ans  =
     
        8.  
     
    -->
    		

    前述の様に、MATLABクローンではFORを利用する事はあまり薦められる事 ではない。寧ろ、次の方法で処理する方が美しく、処理も速い。

    -->z=2*x(find(x>3))
     z  =
     
    !   8.    10.    10.    8. !
    -->2*x.*(x>3)
     ans  =
     
    !   0.    0.    0.    8.    10.    10.    8.    0.    0.    0. !
    
    -->z = 0;
    
    -->tmp=find(x>3); 
     
    -->z(tmp)=2*x(tmp)
     z  =
     
    !   0.    0.    0.    8.    10.    10.    8. !
     
    -->z(tmp)=2*x(x>3)      
     
     
    !   0.    0.    0.    8.    10.    10.    8. !
     
    -->z(x>3)=2*x(x>3)        
     z  =
     
    !   0.    0.    0.    8.    10.    10.    8. !
    
    		

    この例では同じ処理を行っているが、最初の例では3以上の元のみを二倍 にした計算であるのに対し、下の3個の例では配列全体の処理となってい る。尚、scilabのパターンマッチングではT,Fを成分とする行列が得ら れる。これは通常の処理ではMATLABやOctaveと同様の1,0に各々対応する。 その為に2*x.*(x>3)の様な計算も問題なく行えるが、TとFは数値の1と 0では無い為、パターンマッチングによる結果を用いた積で、内部的に1,0 の数値に変換されて演算処理が行なわれる。この処理は必要に応じて 使い分けると良い。次に、2*x(x>3)では成分の指定を特に行っていな いが、一番下の例では成分の指定も入れている。この一番下の例を応用す ればより複雑な入れ換え処理も非常に容易に行える。尚、Octaveでも 同じ処理が行えるが、zに代入を行う場合には予めxとサイズを合せておく 必要がある。以下にOctaveでの例を示す。この例では、Scilabで行った 安易な行列の代入が通用せず、左辺の行列を右辺行列のサイズに合せて おく必要がある事に注目されたい。

    octave:1> x=[1:5,5:-1:1]
    x =
    
      1  2  3  4  5  5  4  3  2  1
    
    octave:2> z(x>3)=2*x(x>3) 
    error: invalid index = 0
    error: evaluating assignment expression near line 2, column 7
    octave:2> z=0
    z = 0
    octave:3>  z(x>3)=2*x(x>3) 
    error: invalid index = 0
    error: evaluating assignment expression near line 3, column 8
    octave:3> z=zeros(size(x));
    octave:4> z(x>3)=2*x(x>3) 
    z =
    
       0   0   0   8  10  10   8   0   0   0
    
    octave:5> 
    		

    前述の様に、未定義の場合や適当な値が入っていても、左辺の代入する 側の大きさに適合していないとエラーになる事を示している。 但し、Scilabの場合は上の処理が自然に行える。又、MATLABでは、 行列のサイズを明記していないと、行ベクトルに勝手に変換される事が あるので、代入を実行する配列は予め配列の大きさを確定する様にした 方が良いだろう。

    この処理の一例としてXが0以上の場合は2、Xが0より小の場合に-1を設定 する必要がある場合、次の様に処理すれば無駄なfor文等のループを使わ ないで簡易に済ます事が可能になる;

    (X>=0)*2+(X<0)*(-1);
    		

    又、この様に処理する事も可能である;

    tmp=(x>=0);
    tmp*2+(1-tmp)*(-1)
    		

    この方法の処理の変種を比較したものを以下に示す。尚、この計算 はPentium !!! 600MHzのPCの結果である。

    octave:1> x=rand(100000,1);
    octave:2> t1=time;(x>=0.5)*2+(x<0.5)*(-1);t2=time;t2-t1
    ans = 0.18966
    octave:3> t1=time;tmp=(x>=0.5);tmp*2+(1-tmp)*(-1);t2=time;t2-t1
    ans = 0.13302
    octave:4> t1=time;tmp=(x>=0.5);tmp*2+tmp-1;t2=time;t2-t1
    ans = 0.12432
    octave:5>x;t1=time;for i1=[1:length(x)] if x(i1)>=0.5; x(i1)=x(i1)*2; else x(i1)=-x(i1);end;end;t2=time;t2-t1
    ans = 33.735
    		

    最初の例ではx>=0.5とx<0.5の二種類のパターンマッチングを実行 している。これに対し、二番目の例でx>=0.5のパターンマッチング のみだが、-1の積を余計に実行している。三番目の例では二番目の例と 似ているが、最後の積を予め実行したもので、最後の例はfor文を用いて、 一々同じ処理を繰り返したものである。最速のものと比較して300倍近く もの差が生じている。結局、手間数では同じ様に見えていても行列の積に 関しては、専門のライブラリが用いられている為に、下手にFOR文を 用いるよりは効率的に計算が行えている。この理由からも、安易なFOR文 の利用は出来るだけ避けて、行列の演算に置換えられるものは徹底して 置換えた方が無難である。

    これと同じ処理をScilabで実行した例を以下に示す。

    --> x=rand(100000,1); 
    
    --> t1=timer();(x>=0.5)*2+(x<0.5)*(-1);t2=timer();t2-t1
     ans  =                                     
     
        0.19  
     
    --> t1=timer();tmp=(x>=0.5);tmp*2+(1-tmp)*(-1);t2=timer();t2-t1
     ans  =
     
        0.3  
     
    --> t1=timer();tmp=(x>=0.5);tmp*2+tmp-1;t2=timer();t2-t1       
     ans  =
     
        0.29  
    		

    このScilabの場合、最良の計算がパターンマッチングを二度行なった もので、他は1.5倍程度の処理時間を必要としている。この事に関しては、 Octaveと比較してScilabではパターンマッチングで生じた行列が Booleanの為、型の変更で時間を消費している為と考えられる。

    -->tmp+(1-tmp);
     
    -->tmp+(tmp-1);
     
    -->tmp+tmp-1;  
          !--error     4 
    undefined variable : %b_a_b                  
     
     
    -->2*tmp-1;    
     
    -->tmp/2-1;
          !--error     4 
    undefined variable : %b_r_s                  
     
     
    -->
    		

    このScilabの例ではtmp=(x>=0.5)として計算を行ってたものであるが、 パターンマッチングで生成したベクトルの演算は基本的に積は問題無いが、 和に関しては(tmp-1)の様な差は許容されるが、tmp+tmpは駄目である。 Boolean型の行列は実数との積や差ではじめて実数に変換され、変換さ れる迄は、通常の四則演算が使えない事に注意する必要がある。その為、 この変換に費やす計算量も考慮する必要がある。

    尚、パターンマッチングに関連して、MATLABやOctaveの便利な命令に anyとall命令がある。先ず、anyは行列データに零でない成分があれば1、 零行列の場合は0を返す。これに対し、allは全ての成分が零でない場合 のみ1を返し、他は0となる。この命令を用いる事で余計な処理を行わず に済む。尚、Scilabにこの命令は無いが自分で作ってしまえば良い。 要するに、sum(a)を判定するプログラムを構築する。anyの場合は sum(a)が0以外の場合のみに1、それ以外は0を返し、allの場合は sum(a)==length(a)の場合のみ1、それ以外は0を返すものとする。 以下にOctaveによる例を示す。

    octave:1> a=rand(4,5)-rand(4,5)
    a =
    
      -0.671536   0.539990   0.205556   0.171495   0.276634
      -0.784795   0.585699  -0.274086  -0.448760   0.131415
      -0.072425   0.276092   0.355440  -0.257676   0.357314
       0.356246   0.061500  -0.318618   0.241485  -0.295221
    
    octave:2> if any(a(:,1)>0)
    > lst=find(a(:,1)>0);
    > b=exp(a(lst,1));
    > end;
    octave:3> b
    b = 1.4280
    		

    上記の様にanyを用いる事で、indexに用いるlstが空リストであるかを 心配する必要が無くなる。この様に一致する個所の存在を判定する事に 利用可能だが、全てが一致するかどうかはanyだけで判別出来ない。 この様に全てが1の場合に1を返す命令にallがある。このallの例を 以下に示す。

    octave:11> all(a(:,1)==a(:,3))
    ans = 0
    octave:12> a(:,1)=a(:,3);
    octave:13> all(a(:,1)==a(:,3))
    ans = 1
    octave:14> 
    		

    この様にMATLAB,Scilabによるパターンマッチングを適用した数値処理 は非常に強力で、下手にfor等のloopを用いるよりも速く、プログラムも 簡潔に済ます事が可能になる。

  7. 便利な行列の定義方法

    等間隔のベクトルの設定を行う場合、数値を直接入力する必要は無い。 Scilabでは[初期値:刻み幅:終値]で設定可能である。

    -->[0:0.1:0.5]
     ans  =
     
    !   0.    0.1    0.2    0.3    0.4    0.5 !
    		

    勿論、for文を用いて次の様に設定する事も可能である。

    -->for i=1:5     
    -->bb(1,i)=0.1*i;
    -->end;
     
    -->bb
     bb  =
     
    !   0.1    0.2    0.3    0.4    0.5 !
     
    		

    但し、Scilab及びMATLABの様な言語では下手に行列処理でforループ を用いると処理速度の大幅な低下の原因となるので注意する。 これはMATLABクローンの内部処理の性質によるもので、for文を用いて 処理を行う場合、指定された変数を一々コピーしてから処理を行う為 である。

    又、対角成分に決まった数値を設定する場合はdiag命令を用いる。 diag命令では対角成分に設定する数値をベクトルで与える。尚、 対角成分を数値分ずらす事も可能である。更に与えられた行列の 対角成分を抜き出す事も可能である。又、対角成分が全て1で他が 全て0の行列の生成はeye命令を用いる。

    
    -->diag([1,2])  
     ans  =
     
    !   1.    0. !
    !   0.    2. !
     
    -->diag([1,2],1)
     ans  =
     
    !   0.    1.    0. !
    !   0.    0.    2. !
    !   0.    0.    0. !
     
    -->diag([1,2],2)
     ans  =
     
    !   0.    0.    1.    0. !
    !   0.    0.    0.    2. !
    !   0.    0.    0.    0. !
    !   0.    0.    0.    0. !
    
    -->diag([1,2],-1) 
     ans  =
     
    !   0.    0.    0. !
    !   1.    0.    0. !
    !   0.    2.    0. !
    
    -->A=rand(3,3)
     A  =
     
    !   0.2808599    0.3374143    0.8808449 !
    !   0.0364674    0.2890952    0.0334877 !
    !   0.8203673    0.0012316    0.726183  !
     
    -->diag(A)
     ans  =
     
    !   0.2808599 !
    !   0.2890952 !
    !   0.726183  !
     
    -->eye(3,2)
     ans  =
     
    !   1.    0. !
    !   0.    1. !
    !   0.    0. !
     
    -->eye(2,3)
     ans  =
     
    !   1.    0.    0. !
    !   0.    1.    0. !
     
    
    		

    上記の例で示す様にrand(m,n)でm行n列の乱数行列を生成する。

    -->rand(3,2)
     ans  =
     
    !   0.0855384    0.6311886 !
    !   0.9449118    0.3323201 !
    !   0.6910084    0.565858  !
     
    -->rand(3,3)
     ans  =
     
    !   0.9333189    0.9648601    0.4600306 !
    !   0.7647275    0.5109398    0.0442728 !
    !   0.3087815    0.4734936    0.3781608 !
    		

    ScilabとMATLABの特徴として、零の成分が多い行列、疎行列 (Sparse matrix)を扱う際に行列の零の個所を記録しない形式が 利用可能である。この疎行列を用いた場合、メモリの消費が押さ えられる利点がある。疎行列の形式を利用した場合、非零の成分 のみを(行,列)の形式で表示する。通常の形式に戻す場合はfull命令 を用いる。疎行列の計算でも通常の形式の行列演算が可能である。 疎行列を利用するとより大きな行列を扱う事が可能になるが、大きな 行列に対して安直にA/B等を計算すると膨大な計算時間を必要とする ので注意を要する。尚、演算の方は通常の行列の演算が利用出来、 通常の行列と混在させて計算を行う事も可能である。

    -->v=[1:4]
     v  =
     
    !   1.    2.    3.    4. !
     
    -->A=diag(v)
     A  =
     
    !   1.    0.    0.    0. !
    !   0.    2.    0.    0. !
    !   0.    0.    3.    0. !
    !   0.    0.    0.    4. !
     
    -->sA=sparse(A)
     sA  =
     
    (    4,    4) sparse matrix
     
    (    1,    1)        1. 
    (    2,    2)        2. 
    (    3,    3)        3. 
    (    4,    4)        4. 
     
    -->sA2=sparse([1:4;1:4]',v)
     sA2  =
     
    (    4,    4) sparse matrix
     
    (    1,    1)        1. 
    (    2,    2)        2. 
    (    3,    3)        3. 
    (    4,    4)        4. 
    
    -->fA=full(sA2)
     ans  =
     
    !   1.    0.    0.    0. !
    !   0.    2.    0.    0. !
    !   0.    0.    3.    0. !
    !   0.    0.    0.    4. !
    
    -->fA-2*sA
     ans  =
     
    ! - 1.    0.    0.    0. !
    !   0.  - 2.    0.    0. !
    !   0.    0.  - 3.    0. !
    !   0.    0.    0.  - 4. !
     
    -->sA-sA2*2
     ans  =
     
    (    4,    4) sparse matrix
     
    (    1,    1)      - 1. 
    (    2,    2)      - 2. 
    (    3,    3)      - 3. 
    (    4,    4)      - 4. 
    
    		

    尚、上記の例でsA2=sparse([1:4;1:4]',v)の個所はMATLABの定義方法 とやや異なる。MATLABでは[1:4;1:4]の様に定義せず、[1:4],[1:4]と 行と列を分けて定義する方式を採用しているので、MATLABとのやり取り には注意が必要である。更に、MATLABクローンのOctaveには標準で SparseMatrixを扱える命令を持たないので、MATLABとプログラムを やり取りする場合には注意する。但し、開発版のOctaveには Sparse Matrix関係の命令も含まれている為、今後に期待出来るだろう。

  8. 多項式の扱い

    MATLAB,OctaveやScilabでは多項式を扱う事も可能である。但し、 数式処理の様に直接多項式を扱えるのはScilabのみである。 MATLABとOctaveでは多項式の係数ベクトルの形で扱う事が可能である。

    octave:2> conv([1,2],[1,-2]) 
    ans =
    
       1   0  -4
    
    octave:3> conv([1,2,2,1],[1,-2])
    ans =
    
       1   0  -2  -3  -2
    
    octave:4>  [aa,bb]=deconv([1,3,3,1],[1,1,1])
    aa =
    
      1  2
    
    bb = -1
    
    octave:5> polyval([1,2,3,4],2)
    ans = 26
    		

    この例では、最初にconv([1,2],[1,-2])で(x+2)*(x-2)を計算している。 結果は[1,0,-4]で、これはx^2-4である。この様にMATLAB/Octaveでは n次の多項式はn+1個の要素を持つリストで表現される。更に、deconv の例ではx^3+3*x^2+3*x+1をx^2+x+1で割った商と余りを計算している。 この場合、x^3+3*x^2+3*x+1=(x^2+x+1)*(x+2)-1である事を示している。 この他にも多項式を扱う関数が幾つか存在するが、どれも直接的な 計算を行うものではなく、あくまでも多項式をリストで扱う。 尚、多項式への具体的な値の代入はpolyvalで行う。この例では多項式 x^3+2*x^2+3*x+4に対し、xに2を代入した例を示している。

    これに対してScilabの方が直接的な利用が行えるが、それでも実際の処理 はMathematicaの様な本格的な数式処理と比べて安易に使えるものではない。 最初に多項式の変数を定義する必要があり、その上、多変数の多項式は 扱えない。以下にそのScilabでの数式の扱い例を示す。

    -->s=poly(0,"s");p1 = s^2-1 
     p1  =
     
             2  
      - 1 + s   
    -->t=poly(0,"t");p2 = t^2-1 
     p2  =
     
             2  
      - 1 + t   
    -->p3 = s^3+1
     p3  =
     
             3  
        1 + s   
     
    -->p1+p3
     ans  =
     
         2   3  
        s + s   
     
    -->p1-p2
         !--error     4 
    undefined variable : %p_s_p                  
     
    -->p1/p3
     ans  =
     
        - 1 + s     
        ---------   
                 2  
        1 - s + s   
     
    -->p1*p3
     ans  =
     
             2   3   5  
      - 1 + s - s + s   
    
    		

    この上の例では、多項式の変数として定義された変数のみが利用可能で ある。多項式同士の演算は同じ変数の多項式に限られるが、和"+"、 差"-"、積"*"、商"/"の計算が行える事を上の例で示している。

  9. 安直な関数の定義

    関数の定義はMATLAB、octaveは殆ど同じで、拡張子が.mのファイル (M-File)を構築すれば良い。もしも、このM-Fileがカレントディレクトリ かM-Fileが存在するディレクトリがpathに含まれていれば、自動的に M-Fileの読み込みと実行を行う。尚、scilabはやや異なり、拡張子が sciのファイルを作っても、自動的に読み込まれる訳では無い。

    このscilabでの命令の直接定義ではdeff関数を用いる。例えば、 与えられた二つのベクトルa,bの各元の積を取る関数としてnekoと 言う名前の関数 の定義を行う必要がある場合;

    -->deff('[z]=neko(x,y)','z=x.*y')
     
    -->neko([1,2,3],[4,3,2])
     ans  =
     
    !   4.    6.    6. !
     
    -->
    		

    の様に、始めに'[z]=neko(x,y)'で関数名と関数の出力値と引数を 設定し、次に実際の処理を設定して行く。この例では'z=x.*y'が 実際の処理である。この例では返す値がzのみであるが、複数の行列 を戻す事も可能である。

    -->deff('[z,a]=neko(x,y)','a=length(x);z=x.*y');
     
    -->[aa,bb]=neko([1,2],[3,4])
     bb  =
     
        2.  
     aa  =
     
    !   3.    8. !
     
    		

    この設定はやや面倒であるが、要するに後述のsciファイルを、 関数名の宣言の個所と、関数本体に分けて記述していると言える。 その為、後側の関数本体にはif文等の通常の制御文を含んでいても、 勿論、問題は無い。一例として、ベクトルの長い方を出力する 関数mikeの例を示す。

    -->deff('a=mike(x,y)','if length(x)>length(y) then a=x;else a=y;end')
     
    -->mike([1,2,3],[1,2])
     ans  =
     
    !   1.    2.    3. !
     
    		

    通常の利用で、一々deffを用いるのは面倒である。この場合、拡張子が sciのファイルを記述し、getfで関数を読み込む。尚、MATLABやoctave と異なり、Scilabが必要な関数を自動で読み込まないので注意する。 尚、getfで関数の読み込みを行う際に、同じ名前の関数が既に定義されて いれば、関数の再定義を行ったと云うWarningが出る。又、ファイル名と 内部に記述された関数名が異なる場合も同様のWarningが出る。 例として、次の関数をneko.sciに記述する。

    function [z]=neko(x,y)
    	if length(x)>length(y) then
     		z = x;
    	elseif length(x) < length(y) then
    		z = y;
    	else
    		z = x./y;
    	end;
    		

    この関数nekoでは長さが等しい場合のみ元同士の商を計算し、 それ以外は長い方のベクトルを出力するものである。この関数の 読み込みはgetfで利用する前に行う必要がある。尚、この例では 関数nekoはdeffの例で既に定義済みのものである。ここで注意する 事に関数名とファイル名はちゃんと合わせておく必要がある。 これは後述のsciファイルやMATLABやOctaveのM-Fileでも同じである。

    以下に、このプログラムを新規に読み込んで実行した例を以下に示す;

    -->getf('WebPage/Math/neko.sci')
    Warning :redefining function: neko                    
     
     
    -->neko([1 2 3],[2 3])
     ans  =
     
    !   1.    2.    3. !
     
    -->neko([1 2],[1 2 3])
     ans  =
     
    !   1.    2.    3. !
     
    -->neko([4 5 6],[1 2 3])
     ans  =
     
    !   4.    2.5    2. !
    		

    これがMATLAB/Octaveの場合にはdeff命令は無く、M-Fileと呼ばれる 拡張子が.mのファイルにプログラムを記述する。MATLAB/Octaveでは M-FileがMATLAB/OctaveのPATH上にあれば、呼び出されると自動的に ロードする。M-Fileの頭は'function [返す変数]=関数名(パラメータ)' となっている点は全く同じである。又、文法的にScilabと非常に良く 似ているが、if文ではthenが不要であり、ややScilabよりは簡易な 言語となっている。例えば上の例をM-Fileに記述すると次の様になる;

    function [z]=nekoneko(x,y)
             if length(x)>length(y)
                z=x;
             else
                z=x./y;
             end;
    		

    ここで、if文の違いからscilabをFORTRAN風、MATLABをC風と言う人も 居るが、MATLABの方がScilabよりは少し砕けた文法になっている程度で 内実は殆ど同じである。

    次に、MATLABやOctaveで関数を直接定義したい場合には、Scilabの様な deff命令を用いずに、M-Fileの内容をそのまま記述すれば良い。 これはX等の環境で、安易にファイルをコピーして貼り付けたもので 問題が無い事を意味する。以下に、Octave上で直接関数を定義している 様子を示す。尚、">"が行の先頭に出ているが、これはOctaveが 返しているプロンプトである。

    octave:1> function [z]=nekoneko(x,y)
    > if length(x)==length(y)
    > z=x./y;
    > else if length(x)>length(y)
    > z=x;
    > else
    > z=y;
    > end
    > end
    > end
    octave:2> nekoneko([1:3],[3:-1:1])
    ans =
    
      0.33333  1.00000  3.00000
    
    octave:3> nekoneko([1:3],[3:-1:0])
    ans =
    
      3  2  1  0
    
    octave:4> nekoneko([1:5],[3:-1:0])
    ans =
    
      1  2  3  4  5
    		

    又、Octaveではfunctionの内容を見たい場合にはtype命令を利用する。 実例を以下に示す。

    octave:6> type nekoneko
    nekoneko is a user-defined function:
    
    function z = nekoneko (x, y)
      if length (x) == length (y)
        z = x ./ y;
      else
        if length (x) > length (y)
          z = x;
        else
          z = y;
        endif
      endif
    endfunction
    octave:7> 
    		

    この様に、functionの記述を見たい場合にはtype function名と 入力する。但し、typeで内容を見る事が可能なものはM-Fileや 利用者定義の変数であり、組込関数の方は見る事が出来ない。 この組込かどうかはwhoで調べる事が可能である。

    octave:8> who
    
    *** currently compiled functions:
    
    length    nekoneko
    
    *** local user variables:
    
    aa  bb
    
    octave:9> type aa
    aa is a user-defined variable
    [ 1, 2, 3 ]
    octave:10> 
    		

    この様に、whoを実行して表示される変数や関数に関してはtypeで内容を 見る事が可能である。尚、pi等の数学的定数は組込関数として扱われて いる。その為、whoを実行しても表には出てこない。

  10. unix命令

    MATLAB、Octave、Scilabでは外部命令を実行する命令が存在する。 Scilabではunx命令を用いる。その使用方法は単純にunix('命令') で行う。以下にunix命令を用いた処理の例を示す;

    -->unix('rm aaaaa')
     ans  =
     
        0.  
     
    -->unix('rm aaaaa')
     ans  =
     
        256.  
    		

    この例ではファイルaaaaaの削除を行なっている。零が返ってくる 場合は問題が無かった事を示す。しかし、256が戻る場合はエラー (この例ではファイルが存在しない)が生じている事を示している。

    因にMATLABの場合は"!"を先頭に付けて外部の命令を起動させる。 !を用いる場合、行末に;を追加すると外部命令のオプションとして 判断される事に注意する必要がある。尚、Octaveでは"!"はMatlab で否定を表わす"not"の意味になるので、MATLABでもOctaveでも動作 するプログラムを開発する場合には特に注意する。このOctaveの場合、 system命令とexec命令の二種類が存在する。exec命令の場合には正常 に処理の実行が終了しするとOctave自体を終了させる。system命令の 場合、実行終了後にoctaveに戻るので、通常はこちらを用いると良い だろう。このsystemを用いた安易な例を示す。

    先ず、systemで起動させるプログラムだが、これは外部のプログラムで、 あれば何でも良い。ここでは次のシェルスクリプトとする。このスクリプト はカレントディレクトリ上でファイルサイズのみをファイルx1に出力する ものである。

    #!/bin/sh
    ls -l | awk '{print $5}'>x1
    		

    次に、systemを用いたプログラムmikeを示す。

    function x = mike ()
      system ("tama");
      load ("x1");
      x = sum (x1);
    endfunction
    		

    このmikeでは先程のシェルスクリプトtamaをsystemで起動させ、 その結果ファイルのx1を読み込み、その総和を計算している。 尚、system命令で実行する命令tamaは環境変数PATHで設定された ディレクトリかカレントディレクトリ上にあれば良い。 又、ディレクトリをsystem("/usr/bin/tama")の様に直接指定しても 良い。更に、system命令で実行するプログラムが引数を必要とする 場合、system("tama mike")の様に通常利用する命令を単に"で括れば 良い。

    例えば、次のシェルスクリプトpochiで指定したディレクトリ上のファイル サイズをx1に出力させる。上述のtamaとの違いは二行目に"$1"が追加され ている点である。

    #!/bin/sh
    ls -l $1 |  awk '{print $5}'>x1
    		

    このpochiに対し、mikeを改良したkuroを以下に示す。

    function x = kuro (wrd)
         evl=["pochi ",wrd];
         system (evl);
         load ("x1");
         x = sum (x1);
    endfunction
    		

    このkuroではevlで文字列の結合を行う。例えば、wrdとして"/usr"が 与えられると、/usr上のファイルサイズの総和を計算する事になるが、 evlでは文字列"pochi "とwrdとして与えられた"/usr"が結合された "pochi /usr"が代入される。この値をsystemで評価(その為にevlを"で 括っていない)して生成されたx1を用いて総和が計算される。

    この様に、systemを用いる事で外部プログラムを用いる事が可能と なり、全てをOctave言語に書き直す手間を省いたり、Octaveへの外部 プログラムの組み込みに関して悩む前に、安易にシステム動作の確認 が行えるので重宝する。

    尚、octaveにはcd,ls及びpwdが利用可能である。働きもUNIXのそれと 全く同じものである。但し、scilabではcdに相当する命令がchdirで chdir('directory')で利用する為、octaveのcdと異なる。しかし、 pwdは同じものだが、lsに相当する命令は無い。

  11. 行列データファイルの読み込み

    Scilab、octave、MATLAB共に行列データの読み込みは比較的容易に 行なえる。先ず、Scilabではread命令で簡易に行列データの読み込み が可能である。このread命令ではファイル名と行列のサイズを設定する。 その為、データを全て読み込む際には行列データの列数を把握する必要 があるが、行数が不明な場合は'-1'を設定すれば良い。又、データの ある列のみ、ある行までのデータの読み込みも可能である。例えば以下 のデータ列を準備する;

           1 2 3
           1 2 3
           1 2 3
           1 2 3
           1 2 3
           1 2 3
           1 2 3
           1 2 3
           1 2 3
           1 2 3
    		

    このデータ列に対し、Scilabのread命令の実行例を以下に示す;

    -->a=read('data',-1,3)
     a  =
     
    !   1.    2.    3. !
    !   1.    2.    3. !
    !   1.    2.    3. !
    !   1.    2.    3. !
    !   1.    2.    3. !
    !   1.    2.    3. !
    !   1.    2.    3. !
    !   1.    2.    3. !
    !   1.    2.    3. !
    !   1.    2.    3. !
     
    -->a=read('data',-1,2)
     a  =
     
    !   1.    2. !
    !   1.    2. !
    !   1.    2. !
    !   1.    2. !
    !   1.    2. !
    !   1.    2. !
    !   1.    2. !
    !   1.    2. !
    !   1.    2. !
    !   1.    2. !
     
    -->a=read('data',-1,1)
     a  =
     
    !   1. !
    !   1. !
    !   1. !
    !   1. !
    !   1. !
    !   1. !
    !   1. !
    !   1. !
    !   1. !
    !   1. !
     
    -->a=read('data',2,3) 
     a  =
     
    !   1.    2.    3. !
    !   1.    2.    3. !
    		

    この例で示す様に-1,3で3列、-1,2で2列のデータ、2,3で先頭の2行を 読み込む。

    MATLABとOctaveではload命令で行列データの読み込み処理が容易に行える。 尚、load命令を利用した場合、データは拡張子を除いたファイル名の変数 に設定される;

    octave:1> load neko.txt
    octave:2> neko
    neko =
    
      1  2  3
      4  5  6
    		

    尚、MATLABとOctaveでload命令を用いる場合、既にファイル名に 対応する変数が存在する場合には-forceオプションを付けないと ファイルの読み込みに失敗する。これはload命令で変数の保護を 行う為である。従って、この様な事態を避ける為、スクリプト ファイルでload命令を記述する場合にはforceオプションを最初 から付けておくと良い。

    octave:1> load neko.txt
    octave:2> neko=1
    neko = 1
    octave:3> load neko.txt
    warning: load: local variable name `neko' exists.
    warning: use `load -force' to overwrite
    error: load: unable to load variable `neko'
    error: evaluating index expression near line 3, column 1
    octave:3> load -force neko.txt
    octave:4> neko
    neko =
    
      1  2  3
      4  5  6
    
    octave:5> 
    		

    この様に単純な行列データの読み込みでは命令一つで容易に行えるが、 MATLAB、Scilab共にはfgets,fscanf等のC風の命令を用いた ファイル操作が可能である。ここで注意する事に、MATLABとOctaveが 非常にC風であるのに対し、Scilabは寧ろFORTRAN風である。 例えば、ファイルのOpen、closeではfopenとfcloseをMATALBとOctaveで 用いるのに対し、Scilabではfileをこれらの処理を主に行う。 又、fscanf等の入力ストリームに対し、実際のファイルの読み書きは readやwriteが主である。更にScilabではfgets命令を持たない。更に、 MATLABとOctaveに限定すると、MATLABのクローンであるOctaveの方が 機能的に優れている。例えば、MATLABでは文字列と数値が混在した ファイルの読み込みにて、fscanfやsscanfを利用する場合はC風 に並びを指定して一気に読み込む事が出来ない。この弱点は、恐ら く出力データが行列に限定されている点によるものだろう。それに対し、 Octaveでは"C"フラグを立てる事で、C風に読み込みが行える。 その際、読み込みの並びに応じて変数を出力する事が可能である。

    次にScilabでは一応、MATLABと似た名前の命令が幾つか用意されて いるが、この命令は名前は似ているが別物と考えた方が良い。 sscanによるデータの変換が追加されているが、C風と言うよりは 基本的にFORTRANのファイル処理の方がより近い。但し、Scilabでは バイナリファイルの書き込みや読み込みも可能であり、テキスト ファイルの処理が主なMATLABやOctaveと比べるとその機能は非常に 高いが、Octave等と比べると作法の違いが目立って多く、分り難い 個所である。

  12. 変数の確認

    MATLABやScilabではシステムが記憶している変数の一覧を見る事が 可能である。この変数の一覧を表示する命令はwho命令である。 例えば、Scilabでwhoを実行すると;

    -->x=1
     x  =
     
        1.  
     
    -->who
    your variables are...
     
     x         startup   ierr      MODE_X    scicos_pal          %helps    
     MSDOS     home      PWD       TMPDIR    percentlib          
     fraclablib          soundlib  xdesslib  utillib   tdcslib   siglib    
     s2flib    roblib    optlib    metalib   elemlib   commlib   polylib   
     autolib   armalib   alglib    intlib    mtlblib   SCI       %F        
     %T        %z        %s        %nan      %inf      old       
     newstacksize        $         %t        %f        %eps      %io       
     %i        %e        
     using       4930 elements  out of    1000000.
              and         46 variables out of       1071
     
     
    -->
    		

    と現在Scilabで定義されている変数が表示される。これがOctaveの 場合は;

    octave:1> x=1
    x = 1
    octave:2> who
    
    *** local user variables:
    
    x
    
    octave:3> 
    		

    となり、ユーザが定義した変数を除いて組込み変数のpi,e,i等は表示 されない。ところが、これらの変数はユーザーが自由に書き換える事の 可能である為、whoで無いからと言って不用意に使わない事が重要である。 尚、whoで表示される変数等の詳細は前述の様にtype命令で確認する事も 可能である。

  13. グラフ表示機能

    グラフ表示機能に関してはGNUPlotを用いるOctaveが明らかに劣って いるが、ScilabとMATLABではそう極端な機能的な差は無い。 実際、ScilabにMATLAB風のグラフを表示させる為のパッケージも存在 している程である。尚、Octaveが利用するGNUPLotは標準では画面を 分割して各々にグラフ表示を行う事が出来ないが、パッチを当てる事 で表示が可能となる様である。但し、Octaveのグラフ表示機能はMATLAB と比べて劣っている事もあって、MATLABの命令がそのまま使えない分野 となっている。とは言え、基本的にグラフの分割表示、曲面のSolid 表示とアニメーションが出来ない程度で、単純な曲線データの表示方法 では両者の使い方に違いは殆ど無い。しかし、Scilabは同様の操作を 行う命令は全く別物である。Scilabのグラフ表示機能は非常に高いが、 引数の設定等が多い。勿論、単純な表示のみであれば容易であり、 Octaveのものと似てなくもない。しかし、その機能を生かそうとすると 覚える事が多く大変である。尚、MATLAB風の命令で同等のグラフ表示 を行う為のScilabパッケージを導入すれば、MATLABとの違いをある程度 吸収する事が可能である。因に、このパッケージはINRIAのScilab関連 のサイトから入手可能である。


数式処理の目次に戻る
目次に戻る