next up previous
Next: Scilabのグラフィックス関数の使い方について Up: 実験方法 Previous: Scilabを終了する

参考プログラム

参考のために, いくつかScilabのプログラムの例を挙げておく。

例 4   オイラー法を用いて2次元数空間における微分方程式

$\displaystyle \frac{d}{dt}\vec{x} = \begin{bmatrix}-0.1 & 1 \ -1 & -0.1 \end{bmatrix} \vec{x}$ (35)

を数値的に解き, 解曲線を2次元空間にプロットすることを考える。 ただし, $ \vec{x}$の初期値を$ [1,1]^T$とする。

周知のように, オイラー法は, 微分方程式(35)を

$\displaystyle \vec{x}_{k+1} = \vec{x}_k+ h \begin{bmatrix}-0.1 & 1 \ -1 & -0.1 \end{bmatrix} \vec{x}_k$ (36)

のように近似する方法である。ここでは$ h=0.1$と取ることにする。

繰り返しの計算には while 文を用いることにする。 また, if文およびbreak文の使い方の例を示すために, やや冗長なのだが,

という制御構造を採用する。さらに, 軌跡の描画のために, 解の履歴をwという 変数に保存しておくことにする。

このためのプログラムは以下の通りである。

A=[-0.1,1;-1,-0.1]; //行列Aの定義
x=[1;1]; //xの初期値の定義
w=x;     //過去の履歴を保存する変数
i=0;     //終了判定用変数
h=0.1;   //ステップ幅
while 1  //無限ループ
        x=x+h*A*x; //オイラー法
        w=[w,x];   //履歴を保存
        i=i+1; //終了判定用変数を増加
        if i>1000 then //終了判定
                break; //ループ脱出
        end //if文終了
end  //while文終了
plot(w(1,:),w(2,:)) //解の軌跡を描画
このプログラムを実行して得られる結果が, 先の図4に挙げたものである。

なお, Scilabのスクリプトにおいて, 文字//およびそれ以降は 注釈とみなされる。注釈は, 人間にとってスクリプトを読みやすく する目的で書かれるもので, プログラム本体とはまったく関係はない。

例 5   関数

$\displaystyle f(x_1,x_2)=x^6+x^2 y^2+y^6$ (37)

の勾配ベクトルは

$\displaystyle \nabla f= \begin{bmatrix}6x^5 + 2x y^2 & 2 x^2 y + 6 y^5 \end{bmatrix}$ (38)

であり, Hesse行列は

$\displaystyle \nabla^2 f= \begin{bmatrix}30 x^4 + 2 y^2 & 4 x y\ 4 x y& 2x^2 + 30y^4 \end{bmatrix}$ (39)

なのであるが, 引数xが与えられたときに これらを求める関数 gradhesseを作ってみよう。 ただし, 関数の定義だけが書かれたファイルfn.sciを 用意することにする。

さて, ファイルfn.sciの内容は以下のようになる。

function g=grad(x)//勾配ベクトル
g(1)=6*x(1)**5 + 2 * x(1)* x(2)**2 
g(2)=6*x(2)**5 + 2 * x(1)**2 * x(2)
g=g'

function h=hesse(x)//ヘッセ行列
h(1,1)=30 * x(1)**4 + 2* x(2)**2
h(1,2)=4  * x(1) * x(2)
h(2,1)=4  * x(1) * x(2)
h(2,2)=30 * x(2)**4 + 2* x(1)**2

ここで定義された関数を使いたいときには, Scilabのウィンドウ内で

getf('fn.sci') $ \hookleftarrow$
と入力する。ファイルがカレントディレクトリにない場合には パスの指定が必要になる。

他のスクリプトから上記で定義された関数 gradと関数hesseを 使いたいときには, スクリプト中のこれらの関数を使用する部分より前のところに

getf('fn.sci') $ \hookleftarrow$
と書いておけばよい。

では, ファイルfn.sciで定義された関数gradを用いて, 2次元空間内で式(38)の勾配ベクトルの正の向き (関数の値が発散する方向)に向かって点を動かしてゆくスクリプトを示そう。

getf('fn.sci');
x=[0.1;-0.1];
w=x;
while(norm(x)<1)
    x=x+grad(x)';
    w=[w x];
end
plot(w)

このスクリプトを実行すると, 図6のようなグラフが得られる。

図 6: 勾配方向に発散してゆく軌道の各座標成分
\includegraphics[width=\textwidth]{divergence.eps}


next up previous
Next: Scilabのグラフィックス関数の使い方について Up: 実験方法 Previous: Scilabを終了する
Shigeru HANBA
平成15年11月16日