MATLABマン
0.目次
【基礎】
1.matlabって何!?
2.文字列を出力(disp文,fprintf文)してみよう!
3.四則演算をしてみよう!
4.変数を使ってみよう!
5.行列を使ってみよう!
5-1.行列の定義
5-2.行列の参照
5-3.行列の四則演算
5-4.多次元行列
6.文字列(string)を扱ってみよう!
6-1.文字列の扱い
6-2.文字列の行列
7.行列をグラフ(plot文)を表示してみよう!
7-1.グラフの基本
7-2.グラフのプロパティ(設定)について
8.条件分岐(if文)を使ってみよう!
8-1.if文って?
8-2.比較演算子について
8-3.else文、elseif文を使おう!
8-4.複数条件の組み合わせ
8-5.論理演算子について
8-6.他の条件分岐(switch文)
9.繰り返し処理(for文)を行ってみよう!
9-1.for文の使い方
9-2.行列とfor文を組み合わせよう!
9-3.FizzBuzz問題
9-4.他の繰り返し(while文)
10.乱数(rand関数、randi関数)で遊ぼう!
10-1.乱数を使ってみよう!(rand関数)
10-2.乱数と行列を組み合わせてみよう!(randi関数)
10-3.乱数を固定してみる!(seed値について)
11.関数(function文)を定義して使ってみよう!
12.変数の保存、読込をやってみよう!(load関数、save関数及びmatファイル)
12-1.save関数
12-2.load関数
13.エラー処理(try文、catch文)でエラーを回避してみよう!
【実践編】
14.運動方程式を計算してグラフに出力してみよう!
15.伝達関数(tf)を使ってみよう!
16.伝達関数の周波数応答をボード線図で表示してみよう!
16-1.その応用
17.状態方程式を組んでみよう!
18.PID制御を使ってみよう!
19.最適制御を使ってみよう!
20.H∞制御を使ってみよう!
【基礎】
1.MATLABって何!?
MATLABとは、制御工学で使う機能をふんだんに詰め込んだプログラミング言語です。
過去に授業で使ったC言語の親戚と考えても差し支えありません。
しかし、MATLABは非常に金額 が高く、個人所有は非常に難しくなっています。
そこで、MATLABを代用できるツールとしてScilabやFreeMatと言ったものがあります。
・MATLAB
MathWorks純正のMATLABです。
お高いですが、その分高性能です。
大きな互換性を持つものにFreeMat等がありますが、
MATLAB同様のソースファイルで完全に代用とできるツールは他にありません。
・Scilab
機能としてはMATLABとほぼ同等のものがあります。
しかし、MATLABに対して互換性がないため、ソースコードの流用は難しくなっています。
・FreeMat
MATLABと強い互換性がありますが、一部の機能に対応できていない、等があります。
しかし、制御を伴わない動作、つまり、本書の基本編までであれば互換性を保つことができます。
以降では、MATLAB、FreeMatを想定した説明を書いていきます。
2.文字列を出力(disp文,fprintf文)してみよう!
まず初めに、文字を表示させてみましょう。
エディタに以下を入力します。
disp('Hello World!')
次にメニューより名前をつけて保存を行います。
この時、大小英数字、アンダーバー(_)以外のものを名前に含めてはいけません。
また、名前の初めに数字がくることもよくありません。
保存したら、これを実行するかF5を押すとコンソールウインドウに、
Hello World!
と表示されていると思います。
また、fprintf文(C言語のprintf文とほぼ同一)を使用することで以下のように書いても同様のものとなります。
fprintf('Hello World!\n')
ここで、\nは改行文字を意味しています。
コンソールで表示されるときに改行として変換されます。
次以降で気が向いたらおまじないとしてエディタの初めに
clc
clear all
close all
と入力しておきましょう。
これを入力すると色々初期化できて散らかることを防いでくれます。
※補足:おまじないのなかみ
clc…コンソールウインドウを初期化
clear all…ワークスペースの初期化
close all…グラフをすべて閉じる
3.四則演算をしてみよう!
四則演算とはいわゆる足し算、引き算、掛け算、割り算のことを言います。
場合によっては累乗や余りの計算を含むこともあります。
ここでは足し算、引き算、掛け算、割り算、累乗及び余りの計算方法を説明します。
足し算、引き算、掛け算、割り算は以下のように計算します。
5+2 %足し算
3-2 %引き算
6*8 %掛け算
6/2 %割り算
2^8 %累乗
これで実行を行うと、以下のように出力され、計算結果が表示されます。
ans =
7
ans =
1
ans =
48
ans =
3
ans =
256
この時の『%』はコメントアウトを意味しており、『%』以降はプログラムに影響を与えない意味のない文字であるということを表しています。
次に、余りについては、10を3で割った余りを計算したいとすると、
rem(10,3)
以上で計算でき、
ans =
1
のようになります。
4.変数を使ってみよう!
次に変数と言うものを使ってみます。
変数とはxとかyみたいな値を代入する文字のことを言います。
代入する文字には関数(disp、fprintf、if、forみたいなの)などで使用されていない文字である必要があります。
また、変数の命名規則は大小英数及びアンダーバー(_)のみ使用可とし、数字が最初にくることは不可とします。
具体的な例とすると、
x8=50;
y0=845;
z1=87426;
このように代入行うことができます。
ここで代入した値は、文字を数字の代わりに入力することで使用することができます。
fprintf('x8=%d\n'x8);
uwaaaaa=x8+y0+z1;
fprintf('uwaaaaa=%d\n',uwaaaaa);
↓
x8=50
uwaaaaa=88321
fprintfの使用方法についてはgoogleさんに尋ねてください。
C言語のprintfと大体同じ使用方法なので、そちらを見ても良いです。(ただし、『"』を使用しないこと)
%dが整数、%fが小数点、%sが文字列です。
また、初めから値がされている変数としてpiがあり、πとして使用することができます。
fprintf('%.6f\n',pi);
↓
3.141593
※補足:『;』について
今まで行末に『;』の記号がついていたことがあったと思います。
MATLABでは、行末に何か付けなければいけないという規則はありませんが、
付けないとコンソール上に変数が表示されてしまうといった仕様があります。
コンソールの上に「出力する」と言う処理は計算に比べても大きく負荷がかかっているため、
表示が必要な場合を除き行末に『;』を付けておく事を推奨します。
5.行列を使ってみよう!
5-1.行列の定義
MATLABでは行列を用いた計算も行うことができます。
行列として代入を行う場合には[]を使用して定義します。
A=[1 2];
これで(1 2)の行列がAに代入されました。
縦に行列を作りたい場合も…
B=[1
2];
または、
B=[1;2];
で作成することができます。
これを組み合わせることで2×2の行列も…
C=[1 2
3 4];
と入力することで定義を行うことができます。
※補足:ちょっと特殊な行列
MATLABの行列には任意の値を代入する意外にも生成する方法があります。
それがzeros、ones、eyeというものです。
・zeros
任意の大きさの零行列を作成します。
disp(zeros(5,8));
↓
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
・ones
全てが1で構成される任意の大きさの行列を作成します。
disp(ones(2,3));
↓
1 1 1
1 1 1
・eye
任意の大きさの単位行列を作成します。
disp(eye(8));
↓
1 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0
0 0 1 0 0 0 0 0
0 0 0 1 0 0 0 0
0 0 0 0 1 0 0 0
0 0 0 0 0 1 0 0
0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 1
5-2.行列の参照
行列を引用したい時は直接代入した文字を用いればよいのですが、
中に入っている一部の値だけを引用したいときがあります。
そういう場合は、以下のように入力します。
例えば、4-1節のCから3を取り出したい場合は…
disp(C(2,1));
↓
3
となります。
C(2,1)は「Cと言う行列変数の2行目の1列目」という意味を表しています。
また、1行目ないしは2行目のみをご所望する場合は『:』を使用することで引用することができます。
1行目の場合だと、
disp(C(1,:));
↓
1 2
となり、C(1,:)は「Cと言う行列変数の1行目の全て」という意味を表します。
列の場合も同様に使用でき、2列目なら、
disp(C(:,2));
↓
2
4
となり、「Cと言う行列変数の全ての2列目」と意味になります。
5-3.行列の指定代入
5-2節で行列の参照方法について説明しました。
これを用いて逆に代入を行うこともできます。
例として元々の行列の一部を書き換えたとすると、
F=[1 2 3 4 5 6];
F(5)=0;
disp(F);
↓
1 2 3 4 0 6
このようにできます。
5-4.行列の一括定義
行列の定義方法で連番を代入する方法があります。
次の二つの方法で行列を定義すると、
A=[1 2 3 4 5 6 7 8 9];
B=1:9;
これらは同じ意味を持ちます。
つまり、『:』演算子は範囲を1刻みの行列として指定することができます。
また、刻み量(ステップ量)を変更したい場合は刻み量を0.1としたとき、
B=1:0.1:9;
と指定することで行えます。
※補足:『:』演算子を用いた行列範囲指定
『:』演算子を使用することで、既存の行列の範囲を指定して参照することができます。
A=1:100;
disp(50:60);
↓
50 51 52 53 54 55 56 57 58 59 60
5-5.行列の四則演算
ここでは行列を用いた計算を行います。
方法としては簡単で、普通の計算と同様に演算子(記号)を使用するだけです。
A=[1 2
3 4];
B=[-4 3
-2 1];
disp('A+B=');
disp(A+B);
disp('A-B=');
disp(A-B);
disp('A*B=');
disp(A*B);
disp('A/B=');
disp(A/B);
↓
A+B=
-3 5
1 5
A-B=
5 -1
5 3
A*B=
-8 5
-20 13
A/B=
2.5000 -5.5000
5.5000 -12.5000
行列の計算なので行数、列数があっていない場合はもちろんエラーとなりますの。
また、掛け算割り算を足し算引き算と同様に各要素で行いたい場合は『.』を使用することで解決できます。
disp('A.*B=');
disp(A.*B);
disp('A./B=');
disp(A./B);
↓
A.*B=
-4 6
-6 4
A./B=
-0.2500 0.6667
-1.5000 4.0000
5-6.多次元行列
一般に行列は二次元のものとされますが、三次元以上も行列として扱うことができます。
A=[1 2
3 4];
B=[5 6
7 8];
AB(1,:,:)=A;
AB(2,:,:)=B;
disp(AB);
↓
(:,:,1) =
1 3
5 7
(:,:,2) =
2 4
6 8
6.文字列(string)を扱ってみよう!
6-1.文字列の扱い
文字列も数字と同様に扱うことができます。
例えば文字列の代入は数字の場合と同様に、
str='StringS';
とすることができます。
このとき、文字は『'』(シングルクォーテーション)で囲む必要があります。
『"』(ダブルクォーテーション)ではエラーとなります。
また、文字列の結合をする場合は、strcatを使用することで行えます。
disp(strcat('aaa','bbb'));
↓
aaabbb
ただし、FreeMatでは対応していないっぽいので注意が必要です。
6-2.文字列の行列
MATLABで文字列を一行の行列として扱うことができます。
A='abcdefg';
disp(A(3));
↓
c
しかし、このままでは複数の文字列を行列として扱うことはできないので、
その場合はcellstr関数を使うことでセル配列といわれる形に変換します。
ただし、各文字列行列は同じ文字列としなければいけないので、足りない文はスペースで埋める必要があります。
strarr=cellstr([
'abcde ',
'123 ',
'aaaaaa']);
disp(strarr);
↓
[abcde]
[123]
[aaaaaa]
また、セル配列から文字列を戻す場合はcharを使用することで行えます。
disp(char(strarr(2)));
↓
123
7.行列をグラフ(plot文)を表示してみよう!
7-1.グラフの基本
ここではMATLABのキモともいえるグラフのプロット機能について説明します。
グラフで次の行列を表示するとします。
X=1:12; %月
Y=[6.3 5.9 10.4 15.0 20.3 23.4 26.8 27.7 23.2 19.1 14.2 6.7]; %平均気温
これを以下のように入力すると、
figure;
plot(X,Y);
↓
図1:月ごとの平均気温(2014)
のように表示することができます。
このときのfigureは「新たな図の定義」、
plotは「プロットする変数及び書式」をそれぞれ行っています。
plot関数には様々なプロパティ(設定)があり、これについては次節で説明をします。
7-2.グラフのプロパティ(設定)について
グラフのプロパティは大まかに以下の種類があります。
figure:新たにグラフを定義する
plot:プロットする変数及び書式を設定する
title:画像タイトルをつける
xlim:グラフの表示範囲(x軸)
ylim:グラフの表示範囲(y軸)
xlabel:x軸に表示する文字列
ylabel:y軸に表示する文字列
grid:グリッド線の表示設定
実際に使用してみると以下のようにできます。
X=1:12; %月
Y=[6.3 5.9 10.4 15.0 20.3 23.4 26.8 27.7 23.2 19.1 14.2 6.7]; %平均気温
figure;
plot(X,Y,'--r','linewidth',2); %XとYでグラフを書く
% 赤の点線 線の太さは 2px
title('月ごとの平均気温(2014)')
xlim([min(X) max(X)]); %Xの最小値から最大値まで
ylim([min(Y) max(Y)]); %Yの最小値から最大値まで
xlabel('月')
ylabel('平均気温 [℃]')
grid on
↓
図2:図の装飾
8.条件分岐(if文)を使ってみよう!
8-1.if文って?
ここからはいわゆる制御構文と言うものを扱っていきます。
今回はその中で最も基本的とされるif文を扱います。
ここで、if文とは?
if文は"if"と言う文字が表すように「もし、○○だったら」を行う構文です。
基本的な構造は、
if 条件式
条件が正しい時にすること
end
となっており、値(num)が偶数であれば2で割る処理を行うとすると、
num=4;
if rem(num,2)==0
%偶数であれば余りは0と等しくなり、割り切れる。rem関数については三章参照。
%また、次節で後述するが、『==』は「等しい」という意味を持つ。
num=num/2;
end
disp(num);
↓
2
のようになります。
また、numに3を代入した時では最終的に3が表示されることが確認できます。
8-2.比較演算子について
前節で条件式に含められていた『==』は比較演算子というものです。
ここでは、条件比較に使用できる比較演算子の種類についてご説明いたします。
次に左値をA、右値をBとした時の例を示します。
A==B :AはBと等しい
A<B :AはBよりも小さい(AはB未満である)
A>B :AはBよりも大きい
A<=B :AはBと同じかそれよりも小さい(AはB以下である)
A>=B :AはBと同じかそれよりも大きい(AはB以上である)
A~=B :AはBと等しくない
8-3.else文、elseif文を使おう!
if文と同時に使える構文にelse文というものがあります。
このelse文というのはif文の「もし、○○であれば～」の条件に、
当てはまらなかった場合に使用でき、「もし、○○であれば～違えば～」
の「違えば～」の位置に相当します。
次にサイコロの目が4以上かどうかをYes/Noで表示する例を示します。
xi=3; %サイコロの目
%もし、サイコロの目が4以上だったら
if 4<=xi
disp('Yes');
%違えば、
else
disp('No');
end
↓
No
また、else文の「違えば～」のあとにif文の「もし○○であれば～」を加えた、
elseif文と言うものがあります。
これを使用することで、2パターンより大きい多パターンの分岐を行うことができます。
例として、サイコロが5以上でTop、サイコロが3以上でMid、それ以外でLowとなる処理の例を以下に示します。
xi=3; %サイコロの目
%もし、サイコロの目が5以上だったら
if 5<=xi
disp('Top');
%違えば、もし、サイコロの目が3以上だったら
elseif 3<=xi
disp('Mid');
%違えば、
else
disp('Low');
end
↓
Mid
8-4.複数条件の組み合わせ
if文で比較を行うときに複数の条件について比較を行いたいことがあります。
例えば、数学で使用される「a<x<b」や、「x<a,b<x」という表記がそれに当たります。
このときのxの範囲を比較するならば、
・a<x<b
if a<x || x<b
%したいこと
end
・x<a,b<x
if x<a || b<x
%すべきこと
end
という形にすることで、行うことができます。
このときの『||』は論理演算子といわれるものの一つで、「○○または××」という意味を持っています。
また、MATLAB(他の言語でも)では、1つ目の例のようなものを以下のように同時に3つ以上の値を繋いで比較することはできないので注意が必要です。
if a<x<b
%これは出来ません
end
8-5.論理演算子について
ここでは論理演算氏というものについて説明します。
論理演算子にはAを一つ目の条件式、Bを二つ目の条件式としたときの例を示します。
A && B :AとBの両方が正しいならばその時の条件式は正しい(論理積)
A || B :AまたはBのどちらかが正しいのならばその時の条件式は正しい(論理和)
また、論理演算子にはもう一つ「否定」という意味を持つものがあり、
次のように入力することで条件式の結果がまったく反対のものとなります。
if ~(1==1)
disp('Yes');
else
disp('No');
end
↓
No
「1==1」の条件は1と1は必ず等しくなるので、正しくなくてはおかしいものです。
しかし、結果としてNoが出てくるということは反対として扱われていることがわかります。
また、これと同様に、
if ~(1==2)
disp('Yes');
else
disp('No');
end
とした場合でも、
Yes
と表示されます。
『~』は条件結果をあべこべにすることが役割ってことですね。
8-6.他の条件分岐(switch文)
ぶっちゃけ、ここはif文でも代用ができるのでこの説については読んでも読まなくても構いません。
条件分岐にはif文だけでなくswitch文というものがあります。
このswitch文とはひとつのものに対して一致するかしないかを判別して処理を行うことのできる構文です。
形は次のようになり、
switch A
case a
Aとaが等しい場合
case b
Aとbが等しい場合
case c
Aとcが等しい場合
otherwise
Aがどのcaseとも等しくない場合
end
このように順次比較を行っていくことができます。
例として、それぞれ変数に入力した値に対して色を表示すると、
color=3;
switch color
case 1
disp('red');
case 2
disp('blue');
case 3
disp('yerrow');
case 4
disp('black');
otherwise
disp('rainbow');
end
↓
yerrow
のようになります。
9.繰り返し処理(for文)を行ってみよう!
9-1.for文の使い方
for文とは繰り返し処理といわれるものを行うための構文です。
この構文は指定した変数の値を順次書き換えることで一定値ごとのカウントを行うことが出来ます。
構文の構成は次のようになります。
for カウンタ変数=初期値:最終値
繰り返す内容
end
例えば、iの変数で1から20までのカウントを行うとすると、
for i=1:20
disp(i);
end
↓
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
また、刻み値を変更したい場合の構文は、
for カウンタ変数=初期値:刻み値:最終値
繰り返したい事柄
end
となり、iを変数としたとき、0から100を5ごとにカウントするときは次のようになります。
for i=0:5:100
disp(i);
end
↓
0
5
10
15
20
25
30
35
40
45
50
55
60
65
70
75
80
85
90
95
100
9-2.行列とfor文を組み合わせよう!
forのカウンタ計算をできる特性を用いて、行列を順次計算させることが可能です。
例として、7章の平均気温を[℃](セ氏温度)から[K](絶対温度)に変換する処理を行います。
C=[6.3 5.9 10.4 15.0 20.3 23.4 26.8 27.7 23.2 19.1 14.2 6.7]; %平均気温
K(1:length(C))=0; %絶対温度
for i=1:length(C)
K(i)=C(i)-273.15;
disp(K(i));
end
↓
-266.8500
-267.2500
-262.7500
-258.1500
-252.8500
-249.7500
-246.3500
-245.4500
-249.9500
-254.0500
-258.9500
-266.4500
※補足:length関数
length(C):Cの行列の個数を取得する。
※補足II:for文を使わない一括計算
今回の例は繰り返しで簡単な計算を行うというものでした。
そのため、今回の例は行列計算を使うことでより簡単に記述することができます。
C=[6.3 5.9 10.4 15.0 20.3 23.4 26.8 27.7 23.2 19.1 14.2 6.7]; %平均気温
K(1:length(C))=0; %絶対温度
K=C-273.15;
disp(K);
MATLABでの繰り返しは他のプログラミング言語に比べても遅くなる傾向があるため、
繰り返しを如何に無くすかということが高速化へのカギとなります。
9-3.FizzBuzz問題
ここで今までの知識を用いてFizzBuzz問題といわれるものを解いてみましょう。
FizzBuzz問題とは1から順番に数字を表示していき、
3の時は「Fizz」、5の時は「Buzz」、3と5の時「Fizz Buzz」とそれぞれ表示するというものです。
これを100まで行ったものを実際に作ってみましょう。
プログラムを作ってみたり、ギブアップしたなら次のボタンを押して確認してみてください。
for i=1:100
if rem(i,3)==0 && rem(i,5)==0
disp('Fizz Buzz');
elseif rem(i,3)==0
disp('Fizz');
elseif rem(i,5)==0
disp('Buzz');
else
disp(i);
end
end
↓
1
2
Fizz
4
Buzz
Fizz
7
8
Fizz
Buzz
11
Fizz
13
14
Fizz Buzz
16
17
Fizz
19
Buzz
Fizz
22
23
Fizz
Buzz
26
Fizz
28
29
Fizz Buzz
31
32
Fizz
34
Buzz
Fizz
37
38
Fizz
Buzz
41
Fizz
43
44
Fizz Buzz
46
47
Fizz
49
Buzz
Fizz
52
53
Fizz
Buzz
56
Fizz
58
59
Fizz Buzz
61
62
Fizz
64
Buzz
Fizz
67
68
Fizz
Buzz
71
Fizz
73
74
Fizz Buzz
76
77
Fizz
79
Buzz
Fizz
82
83
Fizz
Buzz
86
Fizz
88
89
Fizz Buzz
91
92
Fizz
94
Buzz
Fizz
97
98
Fizz
Buzz
9-4.他の繰り返し(while文)
こちらもswitch文と同様に必ずと見る必要は無いので、飛ばしたい方は飛ばしていただいて結構です。
繰り返し構文はfor文の他にwhile文というものがあります。
これは「○○の間～」という意味を持つ構文であり、指定した条件から外れるまでループを行います。
ゆえに、条件が適切でないと無限ループといわれるループから脱出できない現象に陥るため注意が必要です。
構文は次のようになります。
while 条件式
内容
end
たとえば、変数awawaが0から10まで繰り返し処理を行うとすると、
awawa=0
while awawa<=10
disp(awawa);
awawa=awawa+1;
end
↓
1
2
3
4
5
6
7
8
9
10
のようになります。
10.乱数(rand関数)で遊ぼう!
10-1.乱数を使ってみる!
今回は乱数に関して取り扱ってみましょう。
乱数を使うための関数としてはrand関数というものがあります。
これは、0から1を小数点までの乱数を得ることができます。
rand()
↓
ans =
0.8856
これを用いて、0から1までのランダム波形も容易に生成することができます。
rund=rand(1,100);
plot(rund)
↓
図3:ランダム波形
randの()内の値は乱数の生成サイズを表しており、
「1,100」は「1行、100列」の意味となっています。
このとき、()内に一つしか値を入れなかった場合は、
その値が行数、列数両方に適用されます。
10-2.乱数と行列を組み合わせてみよう!
ここでは、整数乱数の生成方法を学びます。
しかし、これはMATLABのバージョン、FreeMatでそれぞれ異なるため、それぞれの方法をここでは説明します。
■MATLAB R2009以降
整数乱数を生成できるrandi関数というものがあります。
これについては当方に確認できる環境がないため割愛します。
下記URLまたは、旧バージョンの方法で適用してください。
http://jp.mathworks.com/help/matlab/ref/randi.html
■FreeMat
上記同様にrandi関数が使用できます。
しかし、MATLABのものとは少々仕様が異なるようでここでは説明しません。
■MATLAB R2008b以前
旧バージョンのMATLABでは、標準で整数乱数を生成できる関数はありません。
そのため、rand関数を変形させることで整数乱数を得ます。
これ以降で例として六面サイコロの処理を作っていきます。
まず、rand関数は0から1までの乱数を生成できるので6倍して範囲をさいころの目にあわせます。
rand()*6
ここで値の切り捨てを行うことで範囲が0から5までの6通りになります。
floor(rand()*6)
これに1を加えることで1から6の範囲の整数乱数が生成できました。
floor(rand()*6)+1
↓
ans =
6
また、整数のみを出力できるという特性より行列と組み合わせて値を選択されることもできます。
omkj=cellstr([ 'daikichi',
'kichi ',
'chukichi',
'shokichi',
'suekichi',
'kyo ',
'daikyo ']);
n=floor(rand()*length(omkj))+1; %length:行列の個数
disp(char(omkj(n))); %文字列の選択&セル配列を文字列に変換
↓
kichi
※補足:値の切り上げ、切捨て、四捨五入
上記で使用したfloor(切捨て)の様に他の小数点処理も行うことができます。
以下にその一覧を示します。
floor :負の方向に対する切捨て
ceil :正の方向に対する切上げ
fix :0の方向に対する切捨て
round :四捨五入
10-3.乱数を固定してみる!(シード値について)
先ほどまでの乱数は実行するたびに値が毎回変化していたはずです。
しかし、シミュレーションを行う場合などでは乱数が毎回一律の方が都合がよいこともあると思います。
ここでは、毎回乱数を一律で生成する方法を説明します。
プログラムにおける乱数は擬似乱数といわれるもので、複雑な計算に基づいて得られています。
その計算に使用されている一つの値がシード(seed)値と言われるもので、
通常ではここに時間に関する値が設定されています。
そのため、シード値は常に変化しているということになります。
しかし、このシード値を"固定化"することで毎回同様の乱数を生成することができます。
■MATLAB R2009以降
rng関数を使用します。
ここでは説明しません。
以下を参照してください。
http://jp.mathworks.com/help/matlab/ref/rng.html
■MATLAB R2008b以前
rand関数をオプション使用することで設定できます。
sdに任意の値・シード値を代入して、以下を記述します。
sd=1; %シード値
rand('seed',sd);
この時のrand関数は乱数の生成は行わず、シード値の値のみを設定しています。
このあとに改めてrand関数を使用することで任意の乱数生成が行えます。
http://jp.mathworks.com/matlabcentral/answers/103815-matlab
■FreeMat
seed関数を使用します。
rand関数を使用する前に以下の記述を行うことで設定できます。
sd=1; %シード値
seed(sd,1);
例:
sd=4649;
rand('seed',sd)
rand()
rand('seed',sd)
rand()
rand()
↓
ans =
0.5136
ans =
0.5136
ans =
0.7890
11.関数(function文)を定義して使ってみよう!
ここでは、関数(function)と言う物を定義してみましょう。
関数とは何かしらの結果を返してくれるプログラムのことを言います。
今まで、私たちが書いてきたプログラムは広義では関数に含まれます。
ここで、この章の題である関数の定義とは、
先ほどまで書いてきたプログラムの1処理をパーツ化することを目的とします。
関数化による利点として、プログラムを小分けすることで整理できたり、
他のプログラムから流用できたり、複雑な処理を見た目上簡易化できたりなど複数あります。
ここから関数の作り方を説明していきます。
ここでの例として10-2節のサイコロを例題としていきます。
まずは、関数にさせたい処理に習って、名前をつけたmファイルを作ります。
サイコロのプログラムなのでsaikoro.mで保存。
このmファイルに次のように記述します。先ほどつけた名前がこの関数の名前となります。
%saikoro.m
function saikoro()
disp('XI[sai]');
end
まずは、関数として呼び出せるかの確認を行います。
関数mファイルと同ディレクトリにmain.mというmファイルを作り、下記を記述し実行します。
%main.m
saikoro()
これで以下のように表示ができれば呼び出し成功です。
XI[sai]
ここで呼び出しを失敗した場合はカレントディレクトリが異なっているなどのディレクトリのトラブルが考えられるので見直してみてください。
呼び出し元と呼び出し先のディレクトリは同一である必要があります。
次にdispを消してサイコロの処理をそのまま記述します。
%saikoro.m
functuon saikoro()
xi=floor(rand()*6)+1;
end
しかし、これでmain.mを動かしても何も表示されないと思います。
ここで戻り値(返り値)というものの設定を行います。
戻り値とは関数を実行したときに返答として返した値のことです。
次のように記述を追加します。
%saikoro.m
functuon xi=saikoro()
xi=floor(rand()*6)+1;
end
↓
4
これでmainに値が帰ってきました。
次に一度に2つのサイコロの値を取得してみます。
%saikoro.m
functuon xi=saikoro()
xi=floor(rand()*6)+1;
sai=floor(rand()*6)+1;
end
しかし、戻り値が一つしかありません。
ここでMATLABの特色として複数の値も行列として値を返すことができます。
次のように記述します。
%saikoro.m
functuon [xi sai]=saikoro()
xi=floor(rand()*6)+1;
sai=floor(rand()*6)+1;
end
また、main.mも次のように書き加えます。
%main.m
[a b]=saikoro()
これを実行することでa,bに値が代入されます。
次は任意の個数のサイコロの結果を取得してみます。
これを行うためにはまず、任意の値を関数に伝える"引数"というものを使用します。
引数とは関数の()の中身の値のことで今まででも多く使ってきたことかと思います。
また、その値が文字列であっても引数は引数です。
ここで新たなファイル、saikoros.mを作り、main.mを次のように書き換えます。
%saikoros.m
function saikoros(num)
disp(num);
end
%main.m
saikoros(3)
それぞれ()内に値が加わりました。
これはmain.mから「3」をsaikoros.mに値が送られており、()内のnumに代入されています。
そのため、main.mを実行すると3が表示されることを確認できるかと思います
これが複数になる場合も『,』でつないで今まで使用してきた引数と同様に記述します。
%saikoros.m
function saikoros(num,nan)
disp(num);
disp(nan);
end
%main.m
saikoros(3,4)
以降ではサイコロ3個と仮定して書いていきます。
%saikoros.m
function saikoros(num)
disp(num);
end
%main.m
saikoros(3)
任意の個数とのことなので繰り返しを用いて行列に値を収めてしまいます。
%saikoros.m
function saikoros(num)
for i=1:num %1からnum(3)まで計算。つまりnum(3)回の計算を行う。
xi(i)=floor(rand()*6)+1;
end
end
MATLABでは行列でも戻り値として返すことができます。
よって、
%saikoros.m
function xi=saikoros(num)
for i=1:num %1からnum(3)まで計算。つまりnum(3)回の計算を行う。
xi(i)=floor(rand()*6)+1;
end
end
main.mを実行すると、
ans =
3 1 6
のようになります。
12.変数の保存、読込をやってみよう!(matファイル及びload関数、save関数)
ここではmatファイルについて説明します。
matファイルとは変数の値について保存を行ったファイルです。
これについてはsave関数を使うことで保存及び編集することができます。
12-1.save関数
save関数を使用しての変数の保存は次の構文で行うことができます。
save('ファイル名','変数名');
または、
save ファイル名 変数名
1000個のランダムな値を保存するときは以下の内何れかで保存することができます。(ファイル名:var_save.mat)
rand_arr=rand(1,1000);
save('var_save','rand');
save('var_save.mat','rand');
save var_save rand
save var_save.mat rand
また、複数変数を保存したい場合にも後ろに引数を継ぎ足すことで可能です。
a=1;
b=2;
save var_save2.mat a b
また、既存のsaveファイルに変数をさらに加える場合は、引数に-appendと加えることで行うことができます。
c=3;
save var_save2 c -append
ここで、matファイルの変数を確認するときはwhos関数を使用して、
whos('-file','var_save2');
または、
whos -file var_save2
で以下のように確認することができます。
Name Size Bytes Class
a 1x1 8 double array
b 1x1 8 double array
c 1x1 8 double array
Grand total is 3 elemants using 24 bytes
12-2.load関数
先ほど保存したmatファイルはload関数を用いることで復元し、変数の読み出しを行うことができます。
構文はsave文とほぼ同様となっています。
matファイルに含まれる全ての変数の読み込みを行う場合は、
load('ファイル名');
または、
load ファイル名
で、行うことができます。
例として12-1節のvarsave.matを読み込む場合は、
load('var_save');
で読み込むことができます。
また、一部変数のみの読み込みを行う場合にも、
load('ファイル名','変数名');
または、
load ファイル名 変数名
とすることで可能です。
13.エラー処理(try文、catch文)でエラーを回避してみよう!
原因の分かっているエラーを回避したいけど、if文で回避させるのが面倒ということがあります。
このような場合にはtry文、catch文を使用することで容易に回避することができます。
構文は次のようになります。
try
エラーが起きるかもしれない処理
catch
エラーが起きたときの処理
end
次に例として変数が存在しない時の処理を行います。
try
disp(z);
disp('OK!');
catch
disp('Not found...');
end
↓
Not found...
存在する場合は、
a=1;
try
disp(a);
disp('OK!');
catch
disp('Not found...');
end
↓
1
OK!
となります。
他のエラーが発生する場合も同様に記述することで弾くことができます。
【実践編】
14.運動方程式を計算してグラフに出力してみよう!
ここからは、実際に運動方程式を用いてシミュレーションを行う方法を勉強していきます。
ここでの例は等速直線運動について取り扱います。
物理学でも勉強したかと思う、次の式を使います。
【速度】
v=v0+a*t [m/s]
【距離】
x=x0+v0*t+a*t2 [m]
ただし、
a:加速度
v:速度
v0:初速度
x:距離
x0:初距離
t:経過時間
この時の経過時間を刻み時間dtとしたときの形に変形すると、
【速度】
vi=vi-1+ai*dt [m/s]
【距離】
xi=xi-1+vi-1*t+ai*dt2 [m]
ただし、
ai:今の加速度
vi:今の速度
vi-1:一つ前の速度
xi:今の距離
xi-1:一つ前の距離
dt:刻み時間
これを用いて400mの直線路を加速し続けた自動車のモデルを考えます。
※補足:今回のプログラムの応用例
ここでは1000mを全力疾走して、全力でブレーキを掛けて停止するチキンレースプログラムを載せておきます。
%% 0-1000 チキンレース
clc;
clear all;
close all;
cnv=60^2/10^3; %[m/s]>>[km/h]:単位換算用変数
a_axl=10.0/cnv; %[km/h/s]>>[m/s^2]:アクセル加速度
a_brk=-35.0/cnv; %[km/h/s]>>[m/s^2]:ブレーキ加速度
l=1000.0; %[m]:コース全長
dt=0.01; %[s]:刻み時間
imax=60/dt; %[-]:最大計算回数。60秒を最高とする。
t(1)=0; %[s]:時間
a(1)=0; %[m/s^2]:加速度
v(1)=0; %[m/s]:速度
x(1)=0; %[m]:走行距離
a_kmh(1)=0; %[km/h/s]:加速度(単位換算)
v_kmh(1)=0; %[km/h]:速度(単位換算)
goal_m=0; %[m]:ゴールまでの残り距離
stop_m=0; %[m]:停止までに必要な距離
%改行
printf('\n');
%[シミュレーションの開始]
for i=2:imax
%[加速度の決定]
%もし、「ゴールまでの残り距離」が「停止までに必要な距離」を超えそうだったら、
%ただし、走行距離が1[m]以上の場合に限る。
if goal_m-1<stop_m && 1<x(i-1)
%ブレーキを踏む(ブレーキの加速度)
a(i)=a_brk; %[m/s]
%違えば、
else
%アクセルを踏む(アクセルの加速度)
a(i)=a_axl; %[m/s]
end
%[運動方程式の計算]
%時間
t(i)=t(i-1)+dt; %[s]
%速度
v(i)=v(i-1)+a(i)*dt; %[m/s]=[m/s]+[m/s^2]*[s]
%距離
x(i)=x(i-1)+v(i-1)*dt+a(i)*dt^2; %[m]=[m]+([m/s]*[s]+[m/s^2]*[s]^2)
%ゴールまでの残り距離
goal_m=l-x(i-1); %[m]=[m]-[m]
%停止までに必要な距離
stop_m=(0^2-v(i)^2)/(2*a_brk); %[m]=([m/s]^2-[m/s]^2)/2*[m/s^2]
%加速度の単位換算
a_kmh(i)=a(i)*cnv; %[m/s^2]>>[km/h/s]
%速度の単位換算
v_kmh(i)=v(i)*cnv; %[m/s]>>[km/h]
%もし、加速度が0以上なら、
if 0<a(i)
%"Accel"と表示(pedalの値)
pedal='Accel';
%違えば、
else
%"Brake"と表示(pedalの値)
pedal='Brake';
end
%[データの表示]
fprintf('%s t=%.2f[s] a=%.2f[km/h/s] v=%.2f[km/h] x=%.2f[m]\r',pedal,t(i),a_kmh(i),v_kmh(i),x(i))
%もし、速度が0以下であるなら、
if v(i)<=0
%改行
printf('\n');
%シミュレーションを終了する(ループを抜ける)
break;
end
%一時停止
pause(0.01);
end
%[各種グラフの表示]
%新たにグラフを定義する
% figure;
%変数をグラフとして表示する
% plot(x,y)
%また、「,'linewidth',2」と最後に加えるとプロット線の太さが変わる
%「2」は線の太さでありここを変更することで自在に変更できる
% plot(x,y,'linewidth',2)
%タイトルをつける
% title('時間-加速度')
%グラフの表示範囲(x軸)
% xlim([min(x) max(x)])
%グラフの表示範囲(y軸)
% ylim([min(y) max(y)])
%x軸に表示する文字列
% xlabel('xラベル')
%y軸に表示する文字列
% ylabel('yラベル')
%グリッド線の表示
% grid on;
%時間-加速度
figure;
plot(t,a_kmh,'linewidth',2)
title('時間-加速度')
xlim([0 max(t)])
ylim([min(a_kmh)-5 max(a_kmh)+5])
xlabel('時間 [s]')
ylabel('加速度 [km/h/s]')
grid on;
%時間-速度
figure;
plot(t,v_kmh,'linewidth',2)
title('時間-速度')
xlim([0 max(t)])
ylim([min(v_kmh) max(v_kmh)+10])
xlabel('時間 [s]')
ylabel('速度 [km/h]')
grid on;
%距離-時間
figure;
plot(x,t,'linewidth',2)
title('距離-時間')
xlim([0 l])
ylim([min(t) max(t)+5])
xlabel('距離 [m]')
ylabel('時間 [s]')
grid on;
%距離-加速度
figure;
plot(x,a_kmh,'linewidth',2)
title('距離-加速度')
xlim([0 l])
ylim([min(a_kmh)-5 max(a_kmh)+5])
xlabel('距離 [m]')
ylabel('加速度 [km/h/s]')
grid on;
%距離-速度
figure;
plot(x,v_kmh,'linewidth',2)
title('距離-速度')
xlim([0 l])
ylim([min(v_kmh) max(v_kmh)+10])
xlabel('距離 [m]')
ylabel('速度 [km/h]')
grid on;