2006 化学反応ダイナミクス概論

ここでは Unix コンソールからの Octave の起動・終了と, 対話コンソールを使った簡単な計算の例を示す.

Octave の起動・終了
"octave" と入力すると Octave のプロンプトが表示される.
```
% octave
octave:1> _
```
Octave を終了するには "exit" または "quit" とタイプする.
```
octave:122> exit
% _
```
簡単な計算
評価したい式を一般的なプログラミング言語と同様に入力するとその結果が表示される. (下線が入力する部分)
```
octave:1> 1.16 * exp(-0.234) / 15
ans = 0.061199
octave:2> ans * 10
ans = 0.61199
octave:3> k = 2
k = 2
octave:4> sqrt(k)
ans = 1.4142
```
- 直前の計算結果は変数 "ans" で参照できる.
- また自分で定義した変数を使うことができる.
行列・ベクトル
行列・ベクトルは, 行の区切りに ";"(セミコロン), 行内の列区切りに " "(空白) または ","(カンマ) を用い, かぎ括弧 "[ ]" でくくる.
次の例は連立一次方程式 AX = B を解く.
```
octave:5> A = [2, 1; 1, 3]
A =
  2  1
  1  3


octave:6> B = [3; 4]
B =
  3
  4


octave:7> X = A \ B
X =
  1
  1
```
逆行列を計算する inv() や転置行列を表す後置演算子 ' などを用いて多種の演算が可能である. (Octave の行列)
```
octave:8> inv(A)
ans =
   0.60000  -0.20000
  -0.20000   0.40000

octave:9> B'
ans =
  3  4

octave:10> B' * B
ans = 25
octave:11> det(A)
ans = 5
```
以下の例では, 行列の固有値と固有ベクトルを計算する.
```
octave:1> A = [-1, 1; 1, -1]
A =
  -1   1
   1  -1

octave:2> [vec, lambda] = eig(A)
vec =
   0.70711   0.70711
  -0.70711   0.70711

lambda =
  -2   0
   0   0
```
微分方程式
Octave には常微分方程式のソルバとして LSODE が搭載されている. 次の例では, ２つの反応で記述されるモデルの連立微分方程式を解く.
d[NO₂]/dt = −k₁[NO₂]
d[O]/dt = k₁[NO₂] − k₂[O₂][O]
1. 解くべき微分方程式を関数 f としてとして定義する. (Octave の関数定義)
```
octave:12> function xdot = f(x, t)
>   xdot(1) = -0.3 * x(1);
>   xdot(2) =  0.3 * x(1) - 20.4 * 210000 * x(2);
> endfunction
```
2. 初期条件 x0 と出力を行う時間 t を設定.
```
octave:13> x0 = [1; 0];
octave:14> t = linspace(0, 10, 51)';
```
  - linspace 関数は一定間隔の数値の並び (横ベクトル) を生成する.
3. lsode を呼び出して計算を実行.
```
octave:15> xout = lsode("f", x0, t);
```
4. 計算結果を, 第一列に時間 t を追加して出力する.
```
octave:16> format short e;
octave:17> [t, xout]
ans =
   0.00e+00   1.00e+00   0.00e+00
   2.00e-01   9.42e-01   6.59e-08
   4.00e-01   8.87e-01   6.21e-08
     :          :          :
   9.80e+00   5.29e-02   3.70e-09
   1.00e+01   4.98e-02   3.49e-09
```
  - format コマンドは出力書式を設定する.
  - [t, xout] は縦ベクトル t と同じ行数の行列 xout を連結した行列である. (Octave の行列操作)
非線形方程式
Octave には非線形方程式を解くための fsolve 関数が用意されている. 次の例では, 非線形方程式 f(x) = 1×10¹³ x^−1.39 exp(−17921 / x) − 1 = 0 を解く.
1. 解くべき非線形方程式を関数 f としてとして定義する. (Octave の関数定義)
```
octave:1> function y = f(x)
>   y = 1e13 * x^(-1.39) * exp(-17921/x) - 1;
> endfunction
```
2. 解の初期推定とともに fsolve を呼び出して解を得る.
```
octave:2> [x, info] = fsolve("f", 800)
x = 873.34
info = 1
```
on-line ヘルプ
Octave のプロンプトで help コマンドを用いるとコマンド・文・関数などのヘルプを参照できる. 例えば上の format コマンドのヘルプを見るには以下のようにする.
```
octave:18> help format
```
ヘルプがページャ (less) で表示される. (less の man を参照すること. ヘルプから抜けるには 'q'.)

Octave 入門