MATLAB/Scilabで理解する数値計算

このページでは「MATLAB/Scialbで理解する数値計算」の関連情報を提供しています.
(2006年8月20日更新)

檜山先生の講義資料のページ

Intel Mac版のScilabが公開されています.
Scilab Home Page
Scilabのインストール
正誤表 (初版)
正誤表 (第2刷)
cos(x)の計算法の改良
簡単な例(モンテカルロ法による円周率の計算)


【Scilab Home Page】

Scilab Home Page はここです.
Scilabバージョン4.0が公開されています.

【Scilabのインストール】

WindowsとLinuxではインストーラが用意されていますので 比較的簡単にインストールできます. この他,Solaris, HP-HX版が公開されています.

MacOSXでもバージョン4.0が使えます.

Scilabは
ここからダウンロードできます(Windows版,Linux版).

MacOSX版:
MacOSX版はここからダウンロードできます. Mac OS X 10.4 版では zip を解凍するだけでOKです.X Window上で動作します.

【Intel MacでのScilab】
Intel Macで従来のPowerPC版のScilabも動作しますが速度が遅くなります.
MacBook Pro 2.0GHz の場合,1000次元の連立一次方程式を解くと, PowerPC版で約4秒,Intel版で約1秒です.
【cos(x)の計算の改良】

110ページに掲載されている mycos2(x) の計算法では,x が pi から 2*pi の範囲のとき -cos(x - pi) を用いていますが,cos(2*pi - x) を 用いることで手間を減らすことができます.

function y = mycos2(x) のなかで
   x(j2) = x(j2) - pi;
となっているところを
   x(j2) = 2*pi - x(j2);
とします.こうすると最後の1行
   y(j2) = -y(j2);
は必要ありません.

同様に sin(x) のときには,
   x(j2) = x(j2) - 2*pi;
のようにします.

x の範囲を pi/4 までとする場合にも,
   -sin(x - pi/2) --> sin(pi/2 - x)
   -cos(x - pi)   --> cos(2*pi - x)
のように用いる式を変えることで同様のことができます.

【モンテカルロ法による円周率の計算 】

MATLABで次のように入力してみましょう.
n = 10;
x = rand(1,n)
y = rand(1,n)
r = x.^2 + y.^2
p = sum(r <= 1)/n*4
plot(x,y,'x');
hold on;
ezplot('cos(t)','sin(t)',[0 pi/2]);
hold off;
何度か繰り返すと異なる値が得られます.


Scilabではグラフのための命令が異なります.
plot2d(x',y',-2);
ezplot('cos(t)','sin(t)',[0 pi/2]);

nを変えて何度も計算するために 次のような関数を用意します.ファイル名は "monte.m"とします.
function p = monte(n)
   x = rand(1,n);
   y = rand(1,n);
   r = x.^2 + y.^2;
   p = sum(r <= 1)/n*4;

MATLAB上で
>> monte(100)
のようにします.100を1000や10000に変えると結果は少しずつ 円周率の値に近づいていきます.引数に同じ値を与えても 乱数が異なるため実行するたびに結果が異なります.
同じ乱数を生成するには
rand('seed',0)
のようにします.0を他の整数に変えると乱数のパターンが変わります.