Home › Monthly Archives › 8月 2013

(20) quiverでベクトル場を描く

始点(x,y)=(1,2)から、XY方向成分(3,2)のベクトルと、
始点(x,y)=(3,5)から、XY方向成分(-1,-4)のベクトルを quiver で描いてみる。

(1)まずは1本ずつ描く場合

octave:1> quiver(1,2,3,2)
octave:2> hold on
octave:3> quiver(3,5,-1,-4)
octave:4> hold off

https://octave.dogrow.net/wp-content/uploads/2013/08/20130828_02_001.jpg

(2)2本のベクトルを行列でまとめて記述する場合

octave:5> v=[1 2 3 2;3 5 -1 -4]
v =
   1   2   3   2
   3   5  -1  -4

octave:6> quiver(v(:,1),v(:,2),v(:,3),v(:,4))

https://octave.dogrow.net/wp-content/uploads/2013/08/20130828_02_002.jpg
両ベクトルの始点は上記(1)と同じだが、上記(1)とはベクトルの長さが異なる…
なぜか?

これは、quiver 内部でベクトル同士が重なり合わないように線の長さ(=scale)を自動調整してくれているから
quiverの第5パラメーターでこのscaleを指定可能。

例えば、両ベクトルの長さを自動調整した長さの30%(scale=0.3)に縮めてみる。

octave:7> quiver(v(:,1),v(:,2),v(:,3),v(:,4),0.3)

https://octave.dogrow.net/wp-content/uploads/2013/08/20130828_02_003.jpg

自動調整をやめさせたい場合は scale=0 を指定すればよい。

octave:8> quiver(v(:,1),v(:,2),v(:,3),v(:,4),0)

https://octave.dogrow.net/wp-content/uploads/2013/08/20130828_02_001.jpg
これで上記(1)と同じように表示された。

(19) csvread()でCSVファイルのデータをロード

こんなCSVファイルを用意した。

[user@dog-server]$ cat test.csv
10,1
20,6
30,18
40,11
50,14
[user@dog-server$

このCSVファイルを csvread でロードする。

octave:5> M=csvread('test.csv')
M =
   10    1
   20    6
   30   18
   40   11
   50   14

octave:6>

ロードしたデータを折れ線グラフで表示してみる。

octave:7> plot(M(:,1),M(:,2))

http://octave.dogrow.net/wp-content/uploads/2013/08/20130828_001.jpg

折れ線グラフではなく、赤色「r」のマーカー「+」でプロットしてみる。

octave:8> plot(M(:,1),M(:,2),'r+')

http://octave.dogrow.net/wp-content/uploads/2013/08/20130828_002.jpg

今度は棒グラフで表示してみる。

octave:9> bar(M(:,1),M(:,2))

http://octave.dogrow.net/wp-content/uploads/2013/08/20130828_003.jpg

(18) 行列計算で中1レベル連立方程式を解く

2元連立1次方程式をOctaveを使って解いてみる。

まずは話を簡単にするために、以下のような簡単な2式の連立方程式を解いてみる。
y = 2x + 4
y = -2x + 16

Octaveでグラフ表示すると以下のような直線が現れる。

octave:21> x=[-10:10]
x =
  -10   -9   -8   -7   -6   -5   -4   -3   -2   -1    0    1    2    3    4    5    6    7    8    9   10

octave:22> y1 = 2 * x + 4
y1 =
  -16  -14  -12  -10   -8   -6   -4   -2    0    2    4    6    8   10   12   14   16   18   20   22   24

octave:23> y2 = -2 * x + 16
y2 =
   36   34   32   30   28   26   24   22   20   18   16   14   12   10    8    6    4    2    0   -2   -4

octave:24> plot(x,y1,x,y2)
octave:25>

まずは与えられた2式を行列で表せるように各項を並び替える。
y = 2x + 4
y = -2x + 16
 ↓
2x – y = -4
-2x – y = -16

これを以下のような行列式で表し、両辺に逆行列をかけて解を求める。

これをOctaveで書くと以下のようになる。

octave:24> A=[2 -1 ; -2 -1]
A =
   2  -1
  -2  -1

octave:25> B=[-4 ; -16]
B =
   -4
  -16

octave:26> Ainv=inv(A)
Ainv =
   0.25000  -0.25000
  -0.50000  -0.50000

octave:27> xy=Ainv * B
xy =
    3
   10

octave:28>

x=3, y=10 が求まった。

※注意
この方法はいつでも使えるわけではない。
行列Aが逆行列を持たない場合、すなわち両方の式の傾きが同じ場合、2本の直線は交点を持たず、連立方程式として成り立たない。
これをOctaveで確認してみる。

octave:34> A=[1 2;4 8]
A =
   1   2
   4   8

octave:35> Ainv=inv(A)
warning: inverse: matrix singular to machine precision, rcond = 0
Ainv =
   Inf   Inf
   Inf   Inf

octave:36>

逆行列が求まらなかった。
2元1次方程式であればイメージがわきやすいが、3元,4元,…,多元連立方程式になると直感的に逆行列の存在有無を判断できない。

次に以下のような4元1次連立方程式を解いてみる。
a + 2b + 3c + 4d = 8
5a + 4b + 6c + 7d = 3
6a + 3b + 4c + 9d = 9
5a + 7b + c + 5d = 5

octave:37> A=[1 2 3 4;5 4 6 7;6 3 4 9;5 7 1 5]
A =
   1   2   3   4
   5   4   6   7
   6   3   4   9
   5   7   1   5

octave:38> B=[8;3;9;5]
B =
   8
   3
   9
   5

octave:39> i=inv(A)
i =
  -5.5674e-01   2.3404e-01   6.7376e-02  -3.5461e-03
   1.6667e-01  -2.5967e-17  -1.6667e-01   1.6667e-01
  -3.1915e-02   3.1915e-01  -1.8085e-01  -9.5745e-02
   3.2979e-01  -2.9787e-01   2.0213e-01  -1.0638e-02

octave:40> A*i
ans =
   1.00000   0.00000   0.00000  -0.00000
   0.00000   1.00000  -0.00000  -0.00000
   0.00000   0.00000   1.00000  -0.00000
   0.00000   0.00000   0.00000   1.00000

octave:41> i*B
ans =
  -3.16312
   0.66667
  -1.40426
   3.51064

octave:42>

a=-3.16312, b=0.66667, c=-1.40426, d=3.51064 が求まった。

★確認
逆行列と元の行列の積は1(=単位行列)となる。上では「octave:40」の A*i でこれを確認している。

(17) rand()とrandn()で乱数生成

rand()は一様分布の乱数生成、randn()は正規分布に従った乱数生成をしてくれる。

(1) まずはrand()で200×200個の乱数を生成し、ヒストグラムでその分布を見てみる。

octave:1> a=rand(200);          % 200x200個の0.0~1.0までの一様乱数を生成
octave:2> b=a*100; c=floor(b);  % 生成した乱数を0~100の整数に変換
octave:3> hist(c)               % ヒストグラム表示

結果は以下のようになった。確かに0~100の全範囲に渡って均等に乱数が生成されている。

(2) 次に、randn()で200×200個の乱数を生成し、ヒストグラムでその分布を見てみる。

octave:1> a=randn(200);         % 200x200個の0.0~1.0までの正規分布乱数を生成
octave:2> b=a*100; c=floor(b);  % 生成した乱数を0~100の整数に変換
octave:3> hist(c)               % ヒストグラム表示

結果は以下のようになった。正規分布に従って乱数が生成されている。

(16) eye()で単位行列を作る

eye()を使うと任意のサイズの単位行列を作ることができる。

octave:1> eye(2)
ans =

Diagonal Matrix

   1   0
   0   1

octave:2> eye(3)
ans =

Diagonal Matrix

   1   0   0
   0   1   0
   0   0   1

octave:3> eye(5)
ans =

Diagonal Matrix

   1   0   0   0   0
   0   1   0   0   0
   0   0   1   0   0
   0   0   0   1   0
   0   0   0   0   1

octave:4> eye(10)
ans =

Diagonal Matrix

   1   0   0   0   0   0   0   0   0   0
   0   1   0   0   0   0   0   0   0   0
   0   0   1   0   0   0   0   0   0   0
   0   0   0   1   0   0   0   0   0   0
   0   0   0   0   1   0   0   0   0   0
   0   0   0   0   0   1   0   0   0   0
   0   0   0   0   0   0   1   0   0   0
   0   0   0   0   0   0   0   1   0   0
   0   0   0   0   0   0   0   0   1   0
   0   0   0   0   0   0   0   0   0   1

octave:5>

(15) ベクトルの内積

Octaveではベクトルの内積も簡単に求められる。
ベクトルaとベクトルbを以下のように定義する。

octave:20> a=[1 2 3]
a =
   1   2   3
octave:21> b=[2 3 4]
b =
   2   3   4

内積は「a・b」だからと「a*b」と書いてしまうとダメ。

octave:22> a*b
error: operator *: nonconformant arguments (op1 is 1x3, op2 is 1x3)

1行3列の行列の積として認識されてしまうのでエラーが出てしまう。

そこで、行列の積が成り立つように形を変えてやる。
すなわち、aが1行3列であれば、bを3行1列に転置してやれば行列の積が成り立つ。
Octaveでは行列の転置を「’」で記述する。

octave:22> b'
ans =
   2
   3
   4
octave:23> a*b'
ans =  20

a=(1,2,3) と b=(2,3,4) の内積は、1×2+2×3+3×4=20 である。

当然だが、以下のように要素の積の総和を求めても同じ。

octave:24> sum(a.*b)
ans =  20

(14) Octave開発バージョンをインストール

Octave開発途中版のインストール手順

■Octaveソースファイルのrepositoryにアクセスするため、Mercurialをインストール。

MercurialのHPはこちら
http://mercurial.selenic.com/wiki/Mercurial

ここからソースファイルをダウンロード
http://mercurial.selenic.com/release/?M=D

HPに書かれた手順でインストール実行
http://mercurial.selenic.com/wiki/UnixInstall

途中、こんなのが出た。

python runrst hgmanpage  --halt warning \
          --strip-elements-with-class htmlonly hg.1.txt hg.1
abort: couldn't generate documentation: docutils module is missing
please install python-docutils or see http://docutils.sourceforge.net/
make[1]: *** [hg.1] Error 255
make[1]: Leaving directory `/home/user/tmp/mercurial-2.7/doc'
make: *** [doc] エラー 2

python-docutilsなるものが必要だそうなのでインストールする。
http://wiki.dreamhost.com/Mercurial

sudo easy_install docutils

■Octave開発版を入手

「Using the Development Sources」に従ってソースファイルを入手
http://www.gnu.org/software/octave/get-involved.html

hg clone http://www.octave.org/hg/octave

■Octave開発版をビルド

手順はココに書かれています。
http://hg.savannah.gnu.org/hgweb/octave/file/tip/etc/HACKING

実行順に以下に列挙します。

cd octave
./bootstrap
mkdir .build
cd .build
../configure
make
make check
$ ./run-octave
$ ./config.status

GUI版Octaveが起動したけれども、期待していたGUIなソースラインデバッグができない…
しかもファイルツリーからファイルを選択してもエディタが起動しない…
機能実装中でリリースはまだまだ先のようです。
使わせていただく立場で好き勝手言って申し訳ありませんが、「開発者の皆様、がんばってください!」
http://octave.dogrow.net/wp-content/uploads/2013/09/gui_octave_20130921.gif