MATLABによるシミュレーション基礎

この資料について

本資料ではシステム制御の学習や研究に必要な動的システムのシミュレーション手法について説明します.シミュレーション環境は基本的にMATLABを対象とします.

シミュレーションとは,ある現象をそれを模擬した装置や計算機を利用して擬似的に再現し,予測や評価に利用することといえます.最近ではデジタル計算機を用いた数値シミュレーションが一般的です. 数値シミュレーションを行うためには,現象の数学的モデルが必要になります.特に工学的な現象の場合,時間に関する微分方程式で表現できることが多いでしょう.その微分方程式を解くことができればシミュレーションが行えます.しかし,微分方程式の解析的な解(数式で表現される解)を求めることは一般的には難しいため,数値解(時間に対する数値の列)を求める手法が用いられます.

この資料ではまず,最も基本的なオイラー法やより精度の高いルンゲ・クッタ法による微分方程式の数値解法をMATLABでのスクリプト例を示しながら説明します.次に,システム制御では基本となる線形システムの状態空間表現について述べ,オイラー法やルンゲ・クッタ法でシミュレ-ションするMATLABスプリプトを示します. 最後に,MATLABのode45関数の利用方法について説明します.

オイラー法によるシミュレーション

一階常微分方程式

微分方程式とは,ある未知関数とその関数の導関数との関係が表されている方程式です.そのうち特に1変数の導関数で表されているものを常微分方程式といいます.システム制御に関わるシミュレーションでは時間に関する導関数で表される場合がほとんどでしょうから,常微分方程式が対象になります.さらに,導関数の階数(次数)が微分方程式の階数となります.

ということで,これから取り扱う最も簡単な微分方程式は以下のような一階の常微分方程式となります.ここで$x(t)$が未知関数,$x_0$はその初期値です.

$$\frac{dx(t)}{dt} = f(x(t), t) \ , \quad x(0) = x_0$$

もちろん実際の問題はより高階の常微分方程式が必要になりますが,この一階の微分方程式の数値解法を応用することで対応できます.

オイラー法

微分方程式の数値解法の基本は,ある時刻$t_k$から微少な刻み幅$\Delta t$(サンプリング時間またはサンプリング周期ともいう)だけ時間がたった後の関数値$x(t_k +\Delta t)$を逐次求めていくことです.これをいかに精度よく(できるだけ少ない計算量で)求めるかが数値解法の善し悪しと言えます.

オイラー法は以下のように$x(t_k +\Delta t)$のテイラー展開の最初の2項だけに基づきます.

$$ \begin{align*} x(t_k + \Delta t) &= x(t_k) + \frac{dx}{dt}\bigg|_{t=t_k} \Delta t + \frac{1}{2} \frac{d^2x}{dt^2}\bigg|_{t=t_k} \ \Delta t ^2 + \cdots \\ &\approx x(t_k) + \frac{dx}{dt}\bigg|_{t=t_k} \Delta t \\ &\approx x(t_k) + f(x(t_k), t_k) \Delta t \end{align*} $$

ここで,表記上の簡便性から,$x(t_k)$を$x_k$と表記することにし,$t_k +\Delta t$を$t_{k+1}$, $x(t_k +\Delta t)$を$x_{k+1}$などと表記することにします.この時オイラー法は以下のように表すことができます.

$$x_{k+1} = x_k + f(x_k, t_k)\Delta t$$

オイラー法による近似をグラフで示すと図のようになります.

オイラー法の概念図

次のサンプリング時刻における$x(t)$の値,すなわち$x_{k+1}$を接線で一次近似して求めていると解釈できます.もとの微分方程式から,接線の傾きが$f(x_k, t_k)$ですから,それに$\Delta t$を乗じた物を$x_k$に加えることで$x_{k+1}$を求めています.図をみると誤差(error)が大きくて,本当に使えるのか心配ですが,$\Delta t$をできるだけ小さくすれば誤差を少なくできます.しかし,$\Delta t$を小さくすると,ある期間(例えば10秒間とか)のシミュレーションをするのにそれだけ計算回数が増えることになるので注意が必要です.

図では$x_k$が真値と一致している場合で書きましたが,実際はそうとも限りません.初期値$x_0$からスタートして,オイラー法を繰り返し使った近似計算は次図のようなイメージになります.$\Delta t$をできるだけ小さくとって,真値から離れすぎないようにする必要があります.なお,この図では近似のための接線を真値のグラフ$x(t)$に接するように引いていますが,$x_k$が真値からずれている場合は図の通りにはなりません.あくまでイメージです.

オイラー法の概念図2

MATLABによる実現

次の微分方程式を例題に,MATLABでオイラー法によるシミュレーションを行うスクリプトを考えてみます.

$$\frac{dx}{dt} = - x + 1, \quad x_0 = 0$$

オイラー法を適用すれば

$$x_{k+1} = x_k +(-x_k + 1)\Delta t$$

とかけるので,繰り返しこの計算を行うMATLABスクリプトの例は以下のようになります.

  • オイラー法による1階微分方程式のシミュレーション
dt = 0.01; % サンプリングタイム  
x_k = 0; % xの初期値  
x = []; % シミュレーション結果格納用変数  
tim =[]; % シミュレーション時刻格納用変数  
  
for t = 0 : dt : 10 % 時刻tは0から10までdt刻み  
    x_k1 = x_k + ( - x_k  + 1) * dt ; % x_{k+1} の計算  
    x = [x x_k]; % xの末尾にx_kを追加  
    tim = [tim t]; % timの末尾にtを追加  
    x_k = x_k1; % 次回の計算のため x_{k} ← x_{k+1}に繰り上げ  
end  
  
plot(tim, x);  
grid on;  
axis([0 10 0 1.5]);  
xlabel('time (s)');  
ylabel('x');

MATLABにおける変数名の都合上,$x_k$をx_k,$x_{k+1}$をx_k1としています. また,シミュレーション結果と時刻を格納する変数をそれぞれxtimとして用意し,x_ktを行の末尾に追加しています.また,forループで繰り返す際に,次回の計算のために今回の$x_{k+1}$を次回の$x_k$として繰り上げておくことが必要です.

高次システムのシミュレーション

高次微分方程式の取り扱い

これまでは1階(次)の微分方程式を扱ってきましたが,現実には2次以上のシステムのシミュレーションが必要になることが多いと思います.その場合は,変数変換をすることで連立多元1階常微分方程式に変形して考えます.

例えば,次のような微分方程式を考えます.これは,長さ$l$の軽い棒の先に質量$m$の重りがついた単振り子を初期角度$\theta_0$の位置から初速度0で放す状況を表しています(下図参照).ここで振り子は空気その他から速度に比例した抵抗を受けるものとし,その係数を$c$としています.

$$\frac{d^2 \theta(t)}{d t^2} = - \frac{g}{l} \sin \theta(t) - \frac{c}{m} \frac{d\theta(t)}{dt} \quad , \theta(0) = \theta_0 \quad , \frac{d\theta(t)}{dt} \bigg|_{t=0} = 0$$

単振り子

このような2次のシステムの場合,次のような2個の変数を新たに定義します.これらを状態変数と呼びます.

$$ \begin{align*} x_1(t) &= \theta(t) \\ x_2(t) &= \frac{d \theta(t)}{dt} \end{align*} $$

状態変数を用いてもとの微分方程式を表すと以下のように1階の微分方程式が2本連立したものになります.

$$ \begin{align*} \dfrac{d x_1(t)}{dt} &= x_2(t) \quad , x_1(0) = \theta_0 \\ \dfrac{d x_2(t)}{dt} &= - \dfrac{g}{l} \sin x_1(t) - \frac{c}{m} x_2(t) \quad , x_2(0) = 0 \end{align*} $$

さらに,2つの状態変数を組にして,

$$ \boldsymbol{x}(t) = \begin{pmatrix}{c}x_1(t)\\x_2(t)\end{pmatrix} \quad , \boldsymbol{x}_0 = \begin{pmatrix}{c}\theta_0\\0\end{pmatrix} $$

のようにベクトル$\boldsymbol{x}(t)$で表現することで,最初の1階常微分方程式と同じ見た目になります.

$$ \frac{d \boldsymbol{x}(t)}{dt} = \boldsymbol{f}(\boldsymbol{x}(t), t) \quad , \boldsymbol{x}(0) = \boldsymbol{x}_0 $$

ここで,$\boldsymbol{f}(\boldsymbol{x}(t), t)$は以下となります.

$$ \boldsymbol{f}(\boldsymbol{x}(t), t) = \begin{pmatrix}x_2(t)\\ - \frac{g}{l} \sin x_1(t) - \frac{c}{m} x_2(t) \end{pmatrix} $$

オイラー法も形式的に次のように表すことができます.

$$ \boldsymbol{x}_{k+1} = \boldsymbol{x}_k + \boldsymbol{f}(\boldsymbol{x}_k, t_k) \Delta t $$

この状態変数表現に基づいたMATLABによるシミュレーション・スクリプトは次のようになります. 基本的には1次のシステムに対するオイラー法のスクリプトと同じですが,x_kx_k1が2行×1列の列ベクトルになっていることが違いです.それに伴い,結果を格納する変数xが最終的に2行×サンプリング点数の行列になります.よって,シミュレーション結果として$\theta$だけをグラフ表示するには1行目だけを取り出す必要があり,そのため,plot関数の引数をx(1,:)としています.

なお,この例ではサンプリングタイム$\Delta t$を0.01秒としていますが,これを0.05秒にするとシミュレーション結果が大きく異なることに気がつくと思います.さらに$\Delta t$を大きくしていくと,それは誤差が大きくなるというようなレベルではなく,全く異なったシミュレーション結果になってしまいます.このようにサンプリングタイムの選択は重要で,基本的には十分小さくとる必要があります.

  • オイラー法による単振り子のシミュレーション
dt = 0.01; % サンプリングタイム  
m = 1; % 質量 1 kg  
l = 0.5; % 振り子長さ 0.5 m  
c = 0.5; % 空気抵抗係数  
x_k = [30/180*pi ;  
         0        ]; % 初期角度30度, 初期速度0  
x = []; % シミュレーション結果格納用変数  
tim =[]; % シミュレーション時刻格納用変数  
  
for t = 0 : dt : 10 % 時刻tは0から10までdt刻み  
    x_k1 = x_k + [               x_k(2)                ;  
                  -9.8/l * sin(x_k(1)) - c/m * x_k(2) ] * dt; % x_{k+1} の計算  
    x = [x x_k]; % xの末尾にx_kを追加  
    tim = [tim t]; % timの末尾にtを追加  
    x_k = x_k1; % 次回の計算のため x_{k} ← x_{k+1}に繰り上げ  
end  
  
plot(tim, x(1, :)); % x_1だけをグラフ表示  
grid on;  
axis([0 10 -1 1]);  
xlabel('time (s)');  
ylabel('\theta (rad)');

ルンゲ・クッタ法によるシミュレーション

2次のルンゲ・クッタ法(ホイン法)

オイラー法は$t_k$における導関数の値$f(x_k, t_k)$を使って$x_{k+1}$を近似していました.ここでは.より精度を高めるために,複数点における導関数の値を利用することを考えます.

先の図のように$t_k$における接線の傾き$f(x_k, t_k)$だけを考えるのではなく,$t_{k+1}$における傾き$f(x_{k+1}, t_{k+1})$も利用して,それらの平均をその区間での傾きとして利用すれば,精度が高まりそうです.その場合,$x_{k+1}$は以下のように考えられます.

$$ x_{k+1} = x_k + \dfrac{1}{2} \left( f(x_k, t_k) + f(x_{k+1}, t_{k+1}) \right) \Delta t $$

しかし,$x_{k+1}$を求める計算に,$x_{k+1}$が含まれているのは無理があるので,$x_{k+1}$をオイラー法で近似した

$$ x_{k+1} = x_k + f(x_k, t_k)\Delta t $$

を代入すれば,

$$ x_{k+1} = x_k + \frac{1}{2} ( f(x_k, t_k) + f(x_k + f(x_k, t_k) \Delta t, t_{k+1}) ) \Delta t $$

とかけます.さらに,逐次計算しやすいように$k_1$,$k_2$という変数を定義して整理すれば,

$$ \begin{align*} k_1 &= f(x_k, t_k) \Delta t \\ k_2 &= f(x_k + k_1, t_{k+1}) \Delta t \\ x_{k+1} &= x_k + \frac{1}{2} ( k_1 + k_2 ) \end{align*} $$

となります.これは2次のルンゲ・クッタ法の一種で,ホイン法または修正オイラー法と呼ばれるものです.2次のルンゲ・クッタ法としては,導関数を取る点として区間の中点を考えた中点法と呼ばれるものもあります.

なお,オイラー法と同様に高階の微分方程式に対しても状態空間表現することでもちろん適用できます.

2次のルンゲ・クッタ法(ホイン法)の原理

前節でオイラー法を用いてシミュレーションした1階の微分方程式をホイン法でシミュレーションするMATLABスクリプトは以下のようになります.微分方程式を複数回計算することになるので,関数fnc()として定義して利用していることに注意してください.

MATLABでユーザ関数を定義するには,functionとendとで定義式を囲み,関数名と同じ名前で拡張子が.mのファイル(この例だとfnc.m)として保存するか,同じスクリプトファイルの最後に記述するかのいずれかの方法をとります.ここでは同じスクリプトファイルに記述しています.

  • ホイン法による1階微分方程式のシミュレーション
dt = 0.01; % サンプリングタイム
x_k = 0; % xの初期値
x = []; % シミュレーション結果格納用変数
tim =[]; % シミュレーション時刻格納用変数

for t = 0 : dt : 10 % 時刻tは0から10までdt刻み
	k1 =  fnc(x_k) * dt; % k1 の計算
	k2 =  fnc(x_k + k1) * dt; % k2 の計算
	x_k1 = x_k + ( k1 + k2) / 2 ; % x_{k+1} の計算
	x = [x x_k]; % xの末尾にx_kを追加
	tim = [tim t]; % timの末尾にtを追加
	x_k = x_k1; % 次回の計算のため x_{k} ← x_{k+1}に繰り上げ
end

plot(tim, x);
grid on;
axis([0 10 0 1.5]);
xlabel('time (s)');
ylabel('x ');

function xdot = fnc(x)  
    xdot = - x + 1 ;  
end

次に2階の微分方程式をホイン法でシミュレーションする例を示します.この場合は$k_1$, $k_2$もベクトルになることに注意してください.

  • ホイン法による単振り子のシミュレーション
dt = 0.01; % サンプリングタイム
x_k = [30/180*pi ;
       0]; % 初期角度30度, 初期速度0
x = []; % シミュレーション結果格納用変数
tim =[]; % シミュレーション時刻格納用変数

for t = 0 : dt : 10 % 時刻tは0から10までdt刻み
	k1 = fnc(x_k) * dt; % k1 の計算
	k2 = fnc(x_k + k1) * dt; % k2 の計算
	x_k1 = x_k + (k1 + k2) /2 ; % x_{k+1} の計算
	x = [x x_k]; % xの末尾にx_kを追加
	tim = [tim t]; % timの末尾にtを追加
	x_k = x_k1; % 次回の計算のため x_{k} ← x_{k+1}に繰り上げ
end

plot(tim, x(1, :)); % x_1だけをグラフ表示
grid on;
axis([0 10 -1 1]);
xlabel('time (s)');
ylabel('\theta (rad)');

function xdot = fnc(x)  
    m = 1; % 質量 1 kg  
    l = 0.5; % 振り子長さ 0.5 m  
    c = 0.5; % 空気抵抗係数  
    xdot = [x(2);
            -9.8/l * sin(x(1)) - c/m * x(2) ] ;  
end

4次のルンゲ・クッタ法

2次のルンゲ・クッタ法では2点における導関数の値を使いました.同様に3点,4点と点数を増やしていくことで原理的には精度がよくなります.しかし実用上,4点を利用した4次のルンゲ・クッタ法で十分なようです.そのため,単にルンゲ・クッタ法というと4次の場合を指すことが多いです.

4次のルンゲ・クッタ法の計算手順は以下のようになります.

$$ \begin{align*} \boldsymbol{k}_1 &= \boldsymbol{f}(\boldsymbol{x}_k, t_k) \Delta t \\ \boldsymbol{k}_2 &= \boldsymbol{f}(\boldsymbol{x}_k + \frac{\boldsymbol{k}_1}{2}, t_k + \frac{\Delta t}{2}) \Delta t \\ \boldsymbol{k}_3 &= \boldsymbol{f}(\boldsymbol{x}_k + \frac{\boldsymbol{k}_2}{2}, t_k + \frac{\Delta t}{2}) \Delta t \\ \boldsymbol{k}_4 &= \boldsymbol{f}(\boldsymbol{x}_k + \boldsymbol{k}_3, t_{k+1}) \Delta t \\ \boldsymbol{x}_{k+1} &= \boldsymbol{x}_k + \frac{1}{6} ( \boldsymbol{k}_1 + 2 \boldsymbol{k}_2 + 2 \boldsymbol{k}_3 + \boldsymbol{k}_4) \end{align*} $$

先の単振り子のシミュレーションをルンゲ・クッタ法で計算するMATLABスクリプトを以下に示します.

  • ルンゲ・クッタ法による単振り子のシミュレーション
dt =  0.01; % サンプリングタイム
m = 1; % 質量 1 kg
l = 0.5; % 振り子長さ 0.5 m
c = 0.5; % 空気抵抗係数
x_k = [30/180*pi ;
0        ]; % 初期角度30度, 初期速度0
x = []; % シミュレーション結果格納用変数
tim =[]; % シミュレーション時刻格納用変数

for t = 0 : dt : 10 % 時刻tは0から10までdt刻み
	k1 = fnc(x_k) * dt; % k1 の計算
	k2 = fnc(x_k + k1/2) * dt; % k2 の計算
	k3 = fnc(x_k + k2/2) * dt; % k3 の計算
	k4 = fnc(x_k + k3) * dt; % k4 の計算
	x_k1 = x_k + (k1 + 2 * k2 + 2 * k3 + k4 ) /6 ; % x_{k+1} の計算
	x = [x x_k]; % xの末尾にx_kを追加
	tim = [tim t]; % timの末尾にtを追加
	x_k = x_k1; % 次回の計算のため x_{k} ← x_{k+1}に繰り上げ
end

plot(tim, x(1, :)); % x_1だけをグラフ表示
grid on;
axis([0 10 -1 1]);
xlabel('time (s)');
ylabel('\theta (rad)');

function xdot = fnc(x)  
    m = 1; % 質量 1 kg  
    l = 0.5; % 振り子長さ 0.5 m  
    c = 0.5; % 空気抵抗係数  
    xdot = [x(2);
            -9.8/l * sin(x(1)) - c/m * x(2) ] ;  
end

MATLAB関数を利用したシミュレーション

ode45関数の使用方法

ルンゲ・クッタ法を用いれば比較でき大きなサンプリングタイムでも精度良くシミュレーションができます.しかし,それでも大き過ぎれば精度は落ちますし,場合によっては発散してしまうこともあります.かといってサンプリングタイムを不必要に小さくすると計算時間が非常に長くなってしまいます.

MATLABにはode45()という関数が用意されています.これは指定した精度に応じてサンプリングタイムを自動的に加減してくれます.しかもサンプリングタイムは一定ではなく,シミュレーション中に必要に応じて小さくなったり大きくなったりします.精度はデフォルトで決まっていますが,オプションとして指定することも可能です.詳しくはode45関数のマニュアルを参照してください.

先の単振り子の自由応答のシミュレーションをode45関数を使って求めるスクリプト例を示します.

  • ode45関数を使った単振り子のシミュレーション1
x0 = [pi/6;
       0];

tspan = [0 10]; % 時間範囲の指定
[t, x] = ode45(@fnc, tspan, x0);

plot(t, x(:, 1));
grid('on');
axis([0 10 -1 1]);
xlabel('time (s)');
ylabel('\theta (rad)');

% 導関数の定義
function xdot = fnc(t, x)
  m = 1;
  l = 0.5;
  c = 0.5;
  xdot = [x(2);
          -9.8/l*sin(x(1))-c/m*x(2)];
end

ode45関数を利用したスクリプトのポイントは以下のようになります.

  • 導関数は時刻tと状態量xをこの順番で引数とする関数として定義する必要があります.$t$を省略することもできません.
  • 定義した導関数の関数名に@マークをつけてode45関数の引数として指定します.
  • シミュレーションの時間範囲を,変数tspanで開始時刻と終了時刻の2要素の行ベクトルとして指定します.
  • サンプリングタイムは自動的に設定されるため,あらかじめ指定する必要はありません.
  • サンプリングタイムは一定ではなく,必要に応じて可変となります.そのため,出力されるシミュレーション結果の状態量xは一定サンプリングではないので注意が必要です.グラフ化する時は同時に出力されるベクトルtとペアにする必要があります.
  • 状態量ベクトルの初期値は列ベクトルで指定しますが,出力される状態量x状態変数が横に並び,時系列に沿って縦長のベクトルになるので注意が必要です.よって特定の状態量だけをプロットするときには列方向に切り出す必要があります.

MATLAB入門

MATLABの利用方法

MATLAB/SIMULINKはMathWorks社が開発・販売する商用の数値解析ソフトウェアである.

東海大学ではMATLAB Campus-Wide Licenseを契約しているので,学生・教職員は個人所有のPCに自由にインストールして利用できる.また,一部のコンピュータ室のPCにはあらかじめインストールされている.さらに,オンライン版ならWebブラウザさえあればどこからでも利用できる(事前にMathWorksアカウントの取得が必要).

東海大学湘南キャンパス湘南校舎情報環境のサイトの,「サポート情報」→「ソフトウェアの利用について」→「MATLABの利用について」に詳細が記載されているので,MathWorksアカウントを取得し,必要に応じて個人PCにインストールするとよい.

その他の数値解析ソフトウェアについて

信号解析・制御系設計のためのCADとして無料で利用できる代表的な数値解析ソフトウェアを以下に示す.

  1. GNU Octave
    MATLABと互換性を持ったフリー(無料)の数値解析ソフトウェアである.MATLABの各種Toolboxに対応するパッケージも有志が開発しており,Octave Packagesとして公開されている.

  2. Scilab フランスのINRIA(国立情報学自動制御研究所)で開発が開始されたフリーの数値解析ソフトウェアである.GNU OctaveほどMATLABとの互換性はないが,SIMULINKに相当する機能も備えている.

  3. Python さまざまな分野で利用出来るオープンソースのプログラミング言語であり,ライブラリ(NumPy, Scipy, Matplotlibなど)を追加することでMATLAB並みの機能が実現出来る.

基本操作

起動

大学コンピュータ室のWindowsPCなら画面下の検索でMATLABを検索して起動する.オンライン版の場合はwebブラウザからMathWorksのサイトにログインして利用する.

カレント・フォルダ

MATLABを起動して開くウィンドウにコマンドを打ち込んでいけば利用はできるが,結果をファイルとして残したり,実験データなとを読み込む場合にはフォルダを意識する必要がある.今現在自分が作業しているフォルダをカレント・フォルダと呼ぶ.

カレント・フォルダを確認するコマンドは,pwdである.cdコマンドでカレント・フォルダを移動できるが,面倒なので,ウィンドウ上部のフォルダーツールバーを利用したほうがいい.なお,オンライン版を利用している場合はクラウドストレージであるMATLAB Drive上にカレント・フォルダがあることになる.必要に応じてPCからファイルをアップロードしたり,PCにダウンロードしたりする.

基本数値・算術コマンド

数値,ベクトル・行列の表現の仕方

MATLAB上でスカラー定数は,

実数:1.6e-31
複素数:1.2+0.8*i (1.2+0.8*jと入力しても良い)

のように入力し,表示もされる.

ベクトル・行列は,

行ベクトル:[1 2 3]
列ベクトル:[1;2;3]
行列:[1 2 3;4 5 6;7 8 9]

のように入力する.要素は実数や複素数,後述の変数(数式)でも構わない.また,要素の区切りはスペース以外にカンマ(,)でもよい.行列は行ベクトルをセミコロン(;)で改行しながら入力するというイメージになる.

変数

MATLAB上では,実数,複素数,行列,ベクトルの区別なく変数に格納することができる.変数名は英文字で始まる英数字の文字列で,デフォルトでは大文字小文字の区別をする.

>> b = 1.2e-3
>> A = [1 2*i;3  1 ]
>> w = -2+3*i

現在宣言されている変数名一覧はコマンドwhoで表示される.

>> who
変数:

A    ans  b    w

定義した変数を削除するには,clearコマンドを使う.引数無しだと,全ての変数をクリアし,

>> clear b w

のように,特定の変数だけをクリアすることもできる.

コマンドライン上での入力と表示

基本的に数値やコマンド・関数をタイプしてenterを押せば,実行される.

>> a = [1 2 3]
a =
  1  2  3

標準状態では結果がすぐ表示されるが,行の最後にセミコロン(;)をつけると表示させないようにすることもできる.

>> a = [1 2 3];

結果が長大なベクトルになる場合が予想されるときなどはセミコロンをつけるとよい.

四則演算,超越関数

四則演算は実数,複素数,行列,ベクトルおよびそれらを内容とする変数に対して実行可能である.ただし,行列の除算A/Bは$S B^{-1}$の意味となる.

sinなどの超越関数も実数,複素数,行列,ベクトルおよびそれらを内容とする変数に対して実行可能であるが,行列,ベクトルに対しては各要素ごとの関数値となる.

>> A = [pi/2  pi;0   pi/6]
>> sin(A)
ans =
  1.00000  0.00000
  0.00000  0.50000

要素ごとの演算

行列やベクトルの和や差はもちろん各要素ごとの和や差となるが,積や除は行列の積の定義に基づいた計算になる.しかし,場合によっては各要素ごとの積や除をとりたい場合がある.その時は,ピリオド(.)を演算子の前につける.

>> A = [1 2;3 4]
>> A*A
ans =
   7  10
  15  22
>> A.*A
ans =
   1   4
   9  16

行列・ベクトルの操作

要素,行,列,小行列の取り出し

行列およびベクトル内の特定の要素を取り出すには,変数名とカッコを用いて次のようにする.

>> a = [1 2 3;4 5 6;7 8 9]
>> a(2,3)
ans = 6

行列の特定の行や列を取り出すには,要素の取り出しの際に列や行番号の代りにコロン(:)を使いる.

>> a = [1 2 3;4 5 6;7 8 9]
>> a(2,:)
ans =
  4     5     6

小行列を取り出すにはコロンを用いて行もしくは列の範囲を指定する.

>> a = [1 2 3;4 5 6;7 8 9]
>> a(2:3,1:2)
ans =
  4     5
  7     8

行列の結合

複数の行列やベクトルをまとめて一つの行列にすることができる.ただし,対応する行数や列数は合わせておく必要がある.

>> a = [1;2;3;4];
>> b = [5;6;7;8];
>> c = [a b]
c =
  1  5
  2  6
  3  7
  4  8

転置行列

行列やベクトルの転置はシングルコーテーション()を用いる.

>> a = [1 2;3 4];
>> a'
ans =
  1     3
  2     4

特定の数列からなるベクトルの作成

MATLABの関数では,シミュレーションの時間軸やボード線図の周波数軸などをベクトルで与えることが多い.そのときは,以下の関数を用いると便利である.

等差数列:
linspace(始値,終値,点数)

たとえば,サンプリングタイム0.01で0から0.99までの点数100のベクトルを作るときは以下のようにする.

>> t = linspace(0,0.99,100);

等比数列:
logspace(log(始値),log(終値),点数)

たとえば,$0.01 = 10^{−2}$から$1000 = 10^3$までを対数分割して100点のベクトルを作るとき.

>> omega = logspace(-2,3,100);

なお,等差数列はlinspace関数を使わなくても次のように:(コロン)だけで作成することもできる.0から1まで0.01刻みで数列を作る場合は以下となる(ただし,ベクトルの点数は101点になることに注意)

>> t = 0:0.01:1;

数値表示のフォーマット

演算結果や変数の内容を表示する桁数を指定できる.標準では5桁の数値で表示される.これでは桁数が少なすぎるとか,指数表示をしてもらいたい時はコマンドformatで指定する.

>> 10/3
ans = 3.3333
>> format long
>> 10/3
ans = 3.33333333333333
>> format longE
>> 10/3
ans =  3.33333333333333e+00

デフォルトの5桁の表示に戻すには,引数無しでformatを実行すればよい.なお,これらはあくまで表示上の変更で,内部計算は8バイト(64 ビット)の倍精度浮動小数点値で行われている.

グラフ

plot関数の基本

1次元のベクトルの組を2次元グラフに表現できる.ベクトルXをx軸,Yをy軸としてグラフ化するには,plot(X, Y)とする.

>> X = linspace(0,10,100);
>> Y = sin(X);
>> plot(X, Y)

線の種類や色を変えたいときは,変数の後に続けてシングルコーテーションでくくったフォーマット文字列を追加する.

>> plot(X, Y, 'b+')

この例は,色をblueとし,実線ではなくて+記号でプロットする.その他の色や線類はhelpで調べられる.

なお,plot関数の引き数としてx軸, y軸, 線の種類を繰り返すことで複数のグラフを重ね書きすることもできる.

>> plot(X1, Y1, 'r-', X2, Y2, 'g-')

その他の装飾

plot関数でグラフを書いた後,xlabel, ylabel, titile関数で軸ラベルとタイトルをつけることができる.

>> xlabel('time (s)')
>> ylabel('sin(x)')
>> title('This is a sample graph')

また,グリッド(格子)をつけるには,

>> grid on

を実行すればよい.

プロットの説明(凡例)はlegend関数で追加できる.特にプロットを複数本重ねる場合は追加すべきである.

 >> legend('sin', 'cos')

のようにプロットした順番に指定する.

不要な場合は,

 >> legend off

で消すことができる.

また,グラフの表示範囲は自動的に設定されるが,複数のグラフを比較するときなど,一律に範囲を指定すべきである.その場合は,axis関数を用いる.

>> axis([xmin xmax ymin ymax])

とすることで,x軸, y軸それぞれの範囲を指定できる.

>> axis auto

で再び自動設定に戻る.

対数グラフ

plot関数のかわりにsemilogx関数を使うと片対数(x軸だけ対数)のグラフになる.両対数はloglog関数を用いる.線種の指定などはplot関数と同じである.

文書への張り込み

グラフをレポートや論文に張り込む場合は,Figureウィンドウの編集メニューから「Figureのコピー」を選べばクリップボードにコピーされるので,それをWordなどのアプリケーション上でペーストすればよい.ただし,オンライン版のMATLABを利用している場合や,印刷に適したベクトル形式の画像フォーマットが欲しい場合は,いったんファイルを経由するほうがよい.その場合はファイルメニューの「保存」で適切なファイル形式を選んで保存する.なお,コマンドウィンドウ上でもprintコマンドを使うことでファイルに出力できる.

>> print('test.png','-dpng')

一つ目の引数はファイル名で,二つ目はファイルの形式を指定するオプションである.-dに引き続いて,pngはPNGフォーマットの意味となる.PNG以外のファイルフォーマットはhelp printで確認できる.

スクリプトファイル

スクリプトファイル(mファイル)とは

これまでのようにコマンドを1行1行打ち込んでいっても作業はできるが,ある程度まとまった処理や繰り返し試行する処理などは,面倒であるし,間違いも見つかりにくい.そこで,その一連の作業をmファイルというテキストファイルとして保存しておくと,実行や修正が非常に楽になる.

mファイルの作成と実行

mファイルはただのテキストファイルなので,MATLAB上のエディタだけでなく,外部の使い慣れたエディタで作成・編集できる.保存する場所は自分のホームフォルダ下に適当な名前のフォルダを作って整理すると良い.

ファイル名は変数名と同じく数字で始まらない英数字とし,拡張子は.mとする.実行はファイル名から.mを抜いた部分をコマンドウィンドウ上で入力すればいい.カレント・フォルダにその名前のファイルがあれば実行される.もしくは,MATLABのツールバーのRunボタンを押すことでも実行される.

path

自作のmファイルをあるフォルダにまとめておき,カレント・フォルダがどこであっても,そのmファイルを実行できるようにすることができる.汎用性の高い関数やツールを作り,将来にわたり使いたい場合に便利である.

path関数によって,現在のコマンド検索pathの表示と変更ができる.

>> path

とすることで現在のpath設定が表示される.これらに,あるフォルダを追加するには以下とする.

>> path(path, '追加したいパス')

この例では,’追加したいパス’で指定したフォルダに自作のmファイルを入れておけば,カレントフォルダがどこであっても実行できるようになる.なお,間違って

>> path('追加したいパス')

とすると,pathが’追加したいパス’だけになり,自作のコマンド以外は全く実行できなくなるので,注意すること.(一度MATLABを終了して再度起動すれば元に戻る)

データファイルの取り扱い

変数の保存と復元

一旦MATLABを終了して次回続きをやりたいときなど,変数の内容をとっておきたいことがある.スクリプトファイルとして作業手順を保存することはできるが,時間のかかる計算が必要な場合,その度に計算し直すのは面倒である. このような時,saveコマンドで変数の内容をファイルに保存できる.

>> save data.mat

とすると,カレントフォルダに全変数の内容を保持したファイルが作成される. このデータを読み込むにはloadコマンドを利用する.

>> clear
>> load data.mat

また,特定の変数を保存するには,以下のように指定する.

>> a = [1 2 3];
>> B = [1 2 3; 4 5 6; 7 8 9];
>> save data.mat a B

しかし,上記の例はMATLABの内部表現(バイナリ形式)であるMATファイル形式での保存であるので,他のアプリケーション(Excelなど)で読み込むことは難しい.そのため,テキスト形式で保存するオプションが用意されている.

>> a = [1 2 3];
>> B = [1 2 3; 4 5 6; 7 8 9];
>> save data.txt a B -ascii

とすれば,ファイルdata.txt に,変数aとBの内容がアスキー形式(テキスト形式)で書き出される.メモ帳など適当なテキストエディタで内容を確認できるはずである.なお,テキスト形式はバイナリ形式にくらべて精度が低く,ファイルサイズも大きくなることに注意が必要である.

タブまたはスペース区切りのテキストデータの読み込み

実験結果をMATLABで処理したりグラフ化する場合,実験用のプログラム(通常はC言語などで自作する)側はMATLABに渡したいデータをテキスト形式のファイルとして保存するのが簡単である.そこで,タブまたはスペース区切りの数列としてテキストファイルで保存しておき,読み込んだMATLAB側で必要に応じて各変数に切り分けることを考える.

例えば,次のような内容のテキストファイルresult.txtがカレントフォルダにあるとする.

1 2 3
2 4 6
3 6 9
4 8 12
5 10 15

このファイルを

>> load result.txt

と読み込むことができる.この時,変数名はファイル名から拡張子をとったものとなる.すなわち,この例の場合resultという名前の変数ができる.もし変数名を指定したい場合は,次のように関数形式で記述する.

>> data = load('result.txt');

先のsaveコマンドで保存したMATファイルの読み込みと異なるのは,一度にひとつの変数しか読み込めないことである.しかし全ての変数を別々のファイルに保存するのは面倒であるので,複数の変数の内容を結合して一つの行列としておき,それを一つのファイルに保存して,あとで切り出すようにするのがよい.

例えば,先の例では,1列目が時間tのデータ,2列目が入力uのデータ,3列目は出力yのデータというように3つの変数を一つの行列にまとめてあると考える.この場合,loadした後で次のように列に切り分ければ,3つの変数として分離できる.

>> data = load('result.txt');
>> t = data(:,1);
>> u = data(:,2);
>> y = data(:,3);

他のアプリケーション向けのテキストデータ書き出し

MATLABでのシミュレーションや解析の結果を他のアプリケーション(例えばExcelなどのスプレッドシートgnuplotといったグラフ化ソフトなど)に渡したい時はテキストファイルとして書き出すのが簡単である.

まず,書き出したい変数を結合して一つの変数にしておく.このときサイズ(行数)が揃っていないとエラーになるので,揃えるか,一緒にすることをあきらめて別のファイルにするか決める必要がある.シミュレーション結果などは,時刻ベクトル,入力ベクトル,出力ベクトルなどが同じ行数のはずなので,まとめて1つの行列にできる.

例えば,時刻t,入力u,出力yの3つの列ベクトルを変数outにまとめてからテキストファイルoutfile.txtに保存するには以下のようにする.

>> out = [t u y];
>> save outfile.txt out -ascii 

なお,テキスト形式での書き出しは,標準では8桁の精度になる.それでは精度が不足する場合は,'-double'オプションを追加すれば,16桁の精度となる.

>> out = [t u y];
>> save outfile.txt out -ascii -double

MATLABスクリプトのいろいろ

キーボード入力,画面への出力

スクリプト実行中にキーボードから数値や文字列を入力させるにはinput命令を用いる.同時にメッセージとして表示させる文字列を指定することができる.

画面への出力はC言語と同じくprintf命令を用いる.

(使用例)

a = input('Please input a: ');
printf('a = %g\n',a);

条件判定

条件判定にはC言語と同じくif命令が利用できる.ただし,C言語のように括弧で条件式は囲まない.また,中括弧({ })は使用せず,endif命令の範囲を指定する.さらに,C言語ではelseifを分けて書くが,MATLABではelseifと続けて書くこと注意.

(使用例)

a = input('Please input a: ');
if a > 0
  printf('a > 0\n');
elseif a > -5
  printf('-5 < a <= 0\n');
else
  printf('a < = -5\n');
end 

繰り返し

while命令

while命令は条件が成り立っている間,命令群を繰り返す.やはり中括弧は用いずに,命令群の終わりはendで指定する.なお,中断はbreak命令である.

while 条件
 命令群
end

(使用例)

sum = 0;
while  sum < 100
 a = input('Please input a: ');
 sum = sum + a;
 printf('sum = %g\n', sum);
end

for命令

for命令はC言語のforを限定的にした感じであり,インデックスとする変数の初期値,増分,終値を指定して命令群を繰り返す.命令群の終わりはend,中断はbreakなのはwhileと同様である.なお,増分を省略すると1と解釈される.

for インデックス変数 = 初期値:増分:終値
 命令群
end

(使用例)

sum = 0;
for  i = 1:2:100
 sum = sum + i;
 if sum >100
  break
 end
end
printf('On i = %d, sum = %g\n', i, sum);

初期値:増分:終値という書き方は,ベクトルの作成の表記と同じである.

ユーザ定義関数

スクリプトファイル内でユーザ独自の関数を定義したり,別ファイルで定義したりできる.

同一ファイル内で定義する場合は以下のように記述する.

function [戻値1, 戻値2, ...] = 関数名(引数1, 引数2, ...)
  % 関数の説明
 命令群
end

MATLABではファイルの末尾に書く必要があり(注:Release 2024a以降ではファイルのどこに書いてもよくなった.ただし,本テキストでは旧バージョンとの互換性を考えてファイルの末尾に書くこととする),Octaveではその関数が呼び出されるより前に記述する必要がある(ただし,ファイルの先頭だと後述のようにファイル名も揃える必要が出てくるので,途中にする必要がある). function命令直後のコメント文はhelpで表示されるので,簡単な説明を入れておくと便利である.

なお,別ファイルで関数を定義する場合は,ファイル名を関数名と同じにした上で,拡張子.mをつける.関数名は既存の関数とかぶらないように,オリジナルなものになるよう注意が必要である.

また,C言語の関数と同様に,関数内で変数はローカル変数となるので,変数名が同じでも内容は別物となることに注意が必要である.

以下の例は,二つのベクトルの各要素ごとに差をとり,それを各々2乗したものの和を計算する関数の例である.

function err = sqerror(a, b)
  % usage: err = sqerror(a, b)
  %  a: column vector
  %  b: column vector, same size as a
  %  err: sum of square errors between a and b
  c  = a - b;
  err = c' * c;
end

TeXTips集

輪講資料などで節番号や数式番号を元文献に合わせたいとき

節番号を合わせる(番号の初期値を変える)

\begin{document}以降に書く.すると,下記の場合は3節から始まる.

\setcounter{section}{2}

数式番号のフォーマットを合わせる

\begin{document}の前に書く.ただし,この場合は節番号を合わせることもしておかないと,数式番号は変わらない.以下は,(節番号.式番号)のようなフォーマットに変更する例である.

\renewcommand{\theequation}{
\thesection.\arabic{equation}}

図の時は、

\renewcommand{\thefigure}{
\thesection.\arabic{figure}}

2段組の文書中で横長の表や図をぶち抜きで入れたい時

表のとき

\begin{table*}[t]
 \begin{tabular}{ccc}
 \end{tabular}
\end{table*}

図のとき

\begin{figure*}[t]
 \begin{center}
  \includegraphics{fig.eps}
  \caption{キャプション}
  \label{}
 \end{center}
\end{figure*}

[b]を指定したい場合は

\usepackage{nidanfloat}

を宣言する.

表のサイズや位置を直したい時

サイズを小さくしたいとき

\begin{table}[ht]
 \begin{center}
  \small                 % ←ここ!
   \begin{tabular}{ccc}
   *****
  \end{tabular}
 \end{center}
\end{table}

上記の位置に\smallなり\scriptsizeを指定する.

位置調整したいとき

\begin{table}[ht]
 \begin{center}
  \hspace*{-2em}              % ←ここ!
  \begin{tabular}{ccc}
  **** 
  \end{tabular}
 \end{center}
\end{table}

例えば,上記のように,上記の場所に\hspace* {**}と書いてあげれば,それに従って表が動きます.ポイントは\hspace*のようにアスタリスクをつけることです.

数式や記号など

ラプラス変換やフーリエ変換の記号を出すには?

筆記体のLを出せばよい.

{\mathcal L}

こんな感じ $\mathcal L$に表示される.逆変換なら

{\mathcal L}^{-1}

とすれば,${\mathcal L}^{-1}$ となる.

フーリエ変換なら

{\mathcal F}

こんな感じ ${\mathcal F}$に表示される.

maketitleで日付を入れたくない場合

maketitleの前に

\date{}

と入れると消える。

デフォルトだと余白が大きくて紙がもったいない.余白を直接指定したい.

geometryパッケージを用いると簡単である. プリアンブルで

\usepackage[top=2cm, bottom=2cm, left=2cm, right=2cm]{geometry}

のように指定すれば良い. なお,pLaTeX2e 新ドキュメントクラス(jsarticleなど)を使う時は,単位をtruecmのようにtrueをつけること. そうしないと,指定より1.2倍大きくなる.

漢字コードがSJISのファイルをUTF-8に変換したい

研究室のMacにはnkfがインストールされているので,

nkf -w --overwrite *.tex

で変換できる. nkfがインストールされていない場合は,homebrewで以下のようにインストールする.

brew install nkf

TeXゼミ基本編

TeXとは...

TeX(テフ,またはテックと読む)は,組版をするためのソフトウェアです.

組版とは,活版印刷の技法で「原稿の指定に従って,順序・字詰め・行数・字間・行間・位置などを正しく組み上げること」を意味しています.つまり,印刷物を作るために文字などを配置してくれるソフトウェアとなります.

なお,実際はTeXをベースに使いやすくしたLaTeXを利用することになります.

TeX環境の準備

PCにインストールする場合

TeX Wikiのサイトを参考にする.

研究室のMacにはTeXLiveをベースにしたMacTeXがインストールされています.一部のWindows機にはW32TeXがインストールされていますが,あまり管理されていないので注意してください.

クラウド版のTeXを利用する場合

PCにインストールしなくても,クラウド上でLaTeXを利用することができるサービスが複数存在します. 無料で利用できるものとしては,株式会社アカリクが提供しているCloud LaTeXや,Overleafがあります.

保存とコンパイル

ソースファイル

TeXはCプログラムと同様に,まず,ソースファイルをエディタで作成しそれをコンパイルして文書を得ます.

ソースファイルは,ファイル名は英数字で,拡張子は.texとします. TeXLiveは使用する文字コードがユニコード(UTF-8)です.

コマンドラインでのコンパイル

TeXのコンパイルに最低限必要なコマンドは以下の通りです.

  • platex ファイル名.tex : dviファイルに変換

  • dvipdfmx ファイル名 : dviファイルをpdfファイルに変換

統合環境

実際には統合環境を使用すると便利なので好きなものを使いましょう.クラウド版のTeXはそのサイトが用意するオリジナルの統合環境で作業することになります.

PCにTeXLiveをインストールした場合,macOS版にはTeXShop,Windows版の場合はTeXworksという統合環境アプリケーションが同胞されます.これら以外にも例えばVisual Studio CodeでTeXを扱えるようにする設定も可能です.

ここではmacOSの場合のTeXShopでの設定を紹介します.まず,TeXShopの「環境設定」の「設定プロファイル」で「pTeX(ptex2pdf)」を選び,エンコーディングを「Unicode (UTF-8)」にします.

過去の資料などでS-JISのコードで書かれたソースはそのままTeXShopで開くと文字化けします.その場合は,以下のコマンドでソースの漢字コードをUTFに変換すればいいです.

nkf -w --overwrite *.tex

基本構成

注意 : S-JISでは,以下の説明で出て来るバックスラッシュ(\)と円マーク(¥)は内部では同じ文字として扱われますが,文字コードをUTF-8とした場合は.バックスラッシュとは別文字の扱いになります.MacのJISキーボードの場合は.optionキーを押しながらキーを押すことでバックスラッシュが入力できます.なお,TexShopのエディタ上ではキーを押せば自動的にバックスラッシュに変換されます.

  1. 自分のホームディレクトリにTeX用のフォルダを作成して下さい.

  2. 適当なエディタで英数字でファイル名を付けて拡張子を .tex にして下さい.これを先ほどのフォルダに保存して下さい.

  3. \documentclass[twocolumn,12pt]{jsarticle} と記入します.[ ]内は段組,文字サイズ,用紙サイズを指定します.もしjsarticleでエラーが出たらjarticleで試してみてください.

  4. \pagestyle{myheadings} と記入します.

  5. \markboth{}{} と記入します. 内はヘッダに出力したい内容を表示します.左のカッコ内は奇数ページ, 右のカッコ内は偶数ページのヘッダです. 例

    {学籍番号 氏名 日付}{学籍番号 氏名 日付}
    
  6. \begin{document} と記入します. 文書の始まりを意味します.

  7. \end{document} と記入します. 文書の終わりを意味します.

以上で文書作成の準備は完了しました.

文書作成

ここからは\begin{document}\end{document}の間に記入していきます.

本文

本文は普通に書き込めばいいですが,ソースファイル上で改行してもそれは無視されます.改行させたい場合は,2個以上の改行(1つ以上の空行)を入れます.または

\\

のようにバックスラッシュ二つでも改行になりますが,インデント(行の頭が一文字下がる)されないので,段落の区切りには使いません.

見出し

報告書やレポートを書くときに区切りとして節を用います.(ex. 「1. TeXとは...」 とか 「5. 文章作成」など)

\section{ }

{ }内には節名を記入します.

\subsection{ }

{ }内には小節名を記入します.

\subsubsection{ }

{ }内には小小節名を記入します.

また,節に番号を付けたくない場合はsectionsubsectionsubsubsectionと{の間に*を付ける事により消すことができます.

\section*{ }

箇条書き

・付き箇条書き

ソース

\begin{itemize}
    \item あいうえお.
    \item かきくけこ.
\end{itemize}

出力結果

・あいうえお.

・かきくけこ.

番号付き箇条書き

ソース

\begin{enumerate}
    \item さしすせそ.
    \item たちつてと.
\end{enumerate}

出力結果

1.さしすせそ.

2.たちつてと.

図のファイル形式

図はPDF形式,JPEGやPNGなどの形式のファイルで挿入できます.ブロック図やイラストなどドロー系の図はPDF形式で書き出して文書に張り込むとクオリティが高いです.

また,写真やスキャンイメージなどのビットマップ系の画像はJPEGやPNG形式のファイルのままでもTeXで扱えますが,画像の量が多いと,処理に時間がかかる場合があります.その場合はこれらの画像ファイルをあらかじめPDF形式に変換しておくと処理が速くなります.ちなみに,Mac上なら,プレビュー.appで開く事ができるファイル形式(PDF, JPEG, PNGなど)ならすべてPDF形式で書き出すことができます.

図の挿入

はじめに \begin{document}より前に

\usepackage[dvipdfmx]{graphicx}

と記入します. そして図を挿入したい場所に

\begin{figure}[htbp]
    \begin{center}
        \includegraphics[width=0.8\columnwidth]{ファイル名.pdf}
    \end{center}
    \caption{図の名前}
    \label{ラベルの名前}
\end{figure}

と記入します.サイズを変更したい時は3行目の数字を小さくしてください.ただし,実際は指定した位置に挿入されるとは限らず,TeXが最も適切な位置を探して挿入します.htbpの意味は次表のとおりで,文字の順番に優先順位が高くなります.

記号 意味
h 記述した位置に出力
t ページ上に出力
b ページ下に出力
p 独立したページに出力

例として3行3列の表の書き方を示します.

ちなみに,下の4行目| |内のl, c, rはそれぞれ,左詰め,真中寄せ,右詰めの意味です.また,l, c,rの数は書きたい列の数で決まって来ます.下の例は3行3列なので,| |の間にl,c,r合わせて3つ用意しました.多すぎる分には問題ありませんが,実際書きたい列より少ないとエラーが出るので注意してください.

また,\hlineは横線を意味します.

\begin{table}[htbp]
    \begin{center}
        \caption{表の名前}
        \begin{tabular}{|l|c|r|}
        \hline
        セル1 & セル2 & セル3 \\ \hline
        セル4 & セル5 & セル6 \\ \hline
        セル7 & セル8 & セル9 \\ 
        \hline
    \end{tabular}
    \label{表のラベル}
    \end{center}
\end{table} 

出力した結果は以下のようになります.

tableの出力例

数式

文書の途中に書く場合

ソース

伝達関数$G(s) = \frac{1}{1+2s}$を定義する.

出力結果

伝達関数$G(s) = \frac{1}{1+2s}$を定義する.

一行の数式

ソース

\begin{equation}
    G(s) = \frac{1}{1+2s} \label{数式名}
\end{equation}

出力結果

$$G(s) = \frac{1}{1+2s}$$

複数行の数式

ソース

\begin{eqnarray}
    G(s) &=& \frac{1}{1+2s} \\
    G(s) &=& \frac{1}{s^2 + s + 1}
\end{eqnarray} 

出力結果

$$ \begin{aligned} G(s) &= \frac{1}{1+2s} \\ G(s) &= \frac{1}{s^2 + s + 1} \end{aligned} $$

このように&= を囲むとそこで揃えてくれます.

なお,このように数式を書くと数式番号がついてしまいますが,もし付けたくない場合は,環境名のequationeqnarrayの最後にアスタリスク(*)をつけるか,付けたくない数式の後ろに以下の命令を付ける事で解消できます.

\nonumber

表,図,数式の番号について

例えば,「◯◯を表△に示す」,「図△に◯◯を示す」,「◯◯の法則より(△)式が求まる」など,表や図,数式の番号を文章中に使う場合,Wordなどの文章作成ソフトでは数字をベタ打していたかもしれませんが,TeXでは次の命令が利用できます.

\ref{}

{}内には表や図を貼付ける際に指定したラベル(\labelの{}の中)を入れます. これを使用することにより,表や図の配置を変えても,自動で表,図などの番号を変えてくれます.なお,ラベルを参照させるときは,TeXのコンパイルを2回繰り返す必要があります.1回目のコンパイルでラベルの対照表を内部で作り,2回目のコンパイルでそれらが使用されます.コンパイルの回数が足りないと,「図??」のように??が挿入されますので注意してください.

使用例

表\ref{◯◯◯}
図\ref{△△△}
式(\ref{×××}),(\ref{×××})式

周波数応答基礎

周波数応答とは

ゲインと位相

安定で線形なシステムに正弦波入力を加えると,定常状態では出力も正弦波となる.出力の正弦波の周波数は入力と同じになるが,振幅の変化や位相差が発生する.

次の図は,伝達関数$\displaystyle \frac{1}{s^2 + 0.5s + 1}$に$\omega= 0.1, 1, 10$ [rad/s]の3種類の正弦波を入力した時の応答をMATLAB(Octave)のlsim関数で計算した結果である.グラフから,過渡状態はおおよそ開始から10秒強で,その後は定常状態となっていることがわかる.これらのグラフから,この伝達関数の場合定常状態において,

  • 周波数が低い場合は入出力にほとんど差がない
  • 周波数が1 [rad/s]の時は入力より出力の振幅の方が大きく,また,入力に対し出力が1/4周期程度遅れている
  • 周波数が高い場合は出力の振幅が非常に小さくなる

ことがわかる.このように周波数に応じて入出力の振幅の関係や時間差(位相ずれ)が変化する.

二次遅れ系の正弦波応答
num = [1];
den = [1 0.5 1];
sys = tf(num, den);
tim = 0:0.01:20;

u1 = sin(0.1*tim);
u2 = sin(1*tim);
u3 = sin(10*tim);

y1 = lsim(sys, u1, tim);
y2 = lsim(sys, u2, tim);
y3 = lsim(sys, u3, tim);

subplot(3,1,1);
plot(tim, u1, tim, y1)
axis([0 20 -2 2]);
xlabel('time (s)');
ylabel('u, y');
title('\omega = 0.1 (rad/s)');
legend('u', 'y');
grid( 'on');

subplot(3,1,2);
plot(tim, u2, tim, y2)
axis([0 20 -2 2]);
xlabel('time (s)');
ylabel('u, y');
title('\omega = 1 (rad/s)');
legend('u', 'y');
grid( 'on');

subplot(3,1,3);
plot(tim, u3, tim, y3)
axis([0 20 -2 2]);
xlabel('time (s)');
ylabel('u, y');
title('\omega = 10 (rad/s)');
legend('u', 'y');
grid( 'on');

そこで入力の正弦波を$u(t)$,定常状態での出力を$y(t)$とし,それぞれ次式で表されるとする.

$$ \begin{array}{l l l} u(t) &=& A_i \sin \omega t \\ y(t) &=& A_o \sin(\omega t + \phi) \end{array} $$

定常状態での波形の例を図に示す.

入出力の波形

このとき,周波数応答は次のゲインと位相の二つの情報の組で表現される.これらが周波数に応じて変化することになる.

  • ゲイン: $g = \displaystyle \frac{A_o}{A_i}$
  • 位相(位相差): $\phi = -\dfrac{\tau}{T} \times 2\pi \, \mathrm{[rad] , or} -\dfrac{\tau}{T}\times 360^{\circ}$

周波数応答関数

システムのゲインと位相は,上記のように実際に正弦波信号を入力して調べなくても,システムの周波数応答関数(周波数伝達関数)から求めることができる.周波数応答関数とは,伝達関数のラプラス演算子$s$を$j\omega$に置き換えたものである.周波数応答関数は角周波数$\omega$を含んだ複素関数であり,その絶対値と偏角がそれぞれゲインと位相に対応する.

$$ \begin{array}{l l l} g &=& |G(j \omega)| \\ \phi &=& \angle G(j \omega) \end{array} $$

なお,ある複素数$z = a + b j$の絶対値と偏角はそれぞれ以下のように求められる.

$$ \begin{array}{l l l} |z| &=& \sqrt{a^2 + b^2} \\ \angle z &=& \tan^{-1} \dfrac{b}{a} \end{array} $$

MATLAB(Octave)での周波数応答の計算

MATLAB(Octave)での複素数の取扱い

MATLAB(Octave)では複素数は以下のように入力し,そのまま実数と同じように扱える.なお,入力する際は,虚数単位はiでもjでもよい.

1 + 2i

例えば,-1の平方根は

> sqrt(-1)
ans =  0 + 1i

と出力される. しかし,変数としてijを使っていると,虚数単位と混同するので,複素数を表す関数complex()を使用した方が間違いがない.また,絶対値は関数abs(),偏角は関数angle()で求められる(単位はラジアン).以下のこれらの使用例を示す.

> z = complex(1, -1)
z =  1 - 1i
> abs(z)
ans =  1.4142
> angle(z)
ans =  -0.78540

MATLAB(Octave)での有理式の取扱い

伝達関数は一般に有理式(分子・分母が多項式)で表される(むだ時間など一部除く).MATLAB(Octave)では多項式を係数を降べきの順に並べたベクトルで表し,有理式は分子・分母それぞれの係数ベクトルのペアで表現する.例えば,伝達関数$\displaystyle G(s) = \frac{s + 2}{s^2 + 3s + 2}$は,次のように二つのベクトルnumdenのペアで表す.

> num = [1 2];
> den = [1 3 2];

また,MATLAB(Octave)では,多項式の変数に値を代入して計算した結果を返してくれる関数polyval()が用意されている.使用例を以下に示す.このように複素数にも対応するので,周波数応答の計算に利用できる.

> f = [1 2 1];
> polyval(f, 1)
ans =  4
> polyval(f, complex(1, 1))
ans =  3 + 4i

周波数応答(ボード線図)のプロット

周波数応答関数のゲインと位相をそれぞれ別のグラフで表したものがボード線図である.横軸は角周波数で対数目盛りとし,ゲインはdBに換算して表示する. ここではMATLAB(Octave)によるボード線図のプロットについて,以下の二つの方法を紹介する.

  • polyval関数を利用して原理に基づいて計算する方法
  • Control System Toolboxのbode関数を利用する方法

polyval関数を利用した方法

以下にpolyval関数を利用してボード線図をプロットするスクリプト例を示す.

変数minpmaxpにそれぞれ周波数軸の下限と上限を10のベキ数で指定する.この例では$10^{-2}=0.01$[rad/s]から$10^2=100$[rad/s]までを指定している.

グラフ表示のウィンドウをsubplot命令で2行1列に分割し,上側にゲイン線図,下側に位相線図をそれぞれ横軸対数表示であるsemilogx関数を使ってプロットしている.図2に実行結果のグラフを示す.

\frac{1}{s^2 + 0.5s+ 1}のBode線図(polyval関数利用)

polyval関数を利用したボード線図の描画スクリプト}

num = [1]; %伝達関数の分子
den = [1 0.5 1];   %伝達関数の分母
minp = -2; %周波数下限のベキ数
maxp = 2;  %周波数上限のベキ数
n = 100;   %周波数分割数
gain =[];
phase = [];
omega = [];

for k = 1 : n
   w = 10^(minp+(maxp-minp)/(n-1)*(k-1));
   z = polyval(num, complex(0, w))/polyval(den, complex(0, w));
   gain(k) = 20*log10(abs(z));
   phase(k) = angle(z)/3.14159*180;
   omega(k) = w;
end

subplot(2,1,1);
semilogx(omega, gain);
xlabel('Frequency (rad/s)');
ylabel('Gain (dB)');
grid('on');

subplot(2,1,2);
semilogx(omega, phase);
xlabel('Frequency (rad/s)');
ylabel('Phase (deg)');
grid('on');

以下のリストは,MATLAB(Octave)固有の要素ごとの演算を利用してforループを用いないスクリプト例である.ここで,logspace関数は等比数列を出力する関数である.

polyval関数を利用したボード線図の描画スクリプト(forループ無しバージョン)

num = [1]; %伝達関数の分子
den = [1 0.5 1];   %伝達関数の分母
minp = -2; %周波数下限のベキ数
maxp = 2;  %周波数上限のベキ数
n = 100;   %周波数分割数

omega = logspace(minp, maxp, n); %周波数ベクトルの生成
z = polyval(num, complex(0, omega))./polyval(den, complex(0, omega));

gain = 20*log10(abs(z));
phase = angle(z)/3.14159*180;

subplot(2,1,1);
semilogx(omega, gain);
xlabel('Frequency (rad/s)');
ylabel('Gain (dB)');
grid( 'on');

subplot(2,1,2);
semilogx(omega, phase);
xlabel('Frequency (rad/s)');
ylabel('Phase (deg)');
grid('on');

bode関数を利用した方法

MATLABのControl System Toolboxに含まれるbode関数を利用すると,簡単にボード線図を描画できる.以下にスクリプト例と実行結果のグラフを示す.なお,この例では変数omegaで周波数ベクトルを指定しているが,省略することもできる.その場合はbode関数が自動的に周波数範囲を設定して描画する.

\frac{1}{s^2 + 0.5s+ 1}のBode線図(bode関数利用)

bode関数を利用したボード線図の描画スクリプト

num = [1]; %伝達関数の分子
den = [1 1 1]; %伝達関数の分母
minp = -2; %周波数下限のベキ数
maxp = 2;  %周波数上限のベキ数
n = 100;   %周波数分割数

omega = logspace(minp, maxp, n); %周波数ベクトルの生成
sys = tf(num, den); %伝達関数からシステムマトリックスへの変換

bode(sys, omega);

先のリストの例のようにbode関数単体でグラフまで描画するが,ゲイン線図だけを描画したい時や,実験結果と重ねてグラフ化したい時などは,次のように,bode関数の左辺に変数をおくと,グラフ表示は行わず周波数応答の計算結果だけが得られる(ゲインは絶対値,位相は度の単位).

[gain phase omega] = bode(sys, omega);

これを利用したスクリプト例を以下のリストに示す.これはボード線図のゲイン線図だけをプロットする例である.実行結果を図に示す.

なお,bode関数が多入出力に対応しているため,MATLABの場合左辺の変数gainphaseは3次元の行列として出力される(対象システム1入力1出力であっても).よって,plot関数で使用する前にこのスクリプト例のようにsqueeze関数で1次元ベクトルに変換する必要があるので注意する.

\frac{1}{s^2 + 0.5s+ 1}のBode線図(bode関数利用,ゲイン線図のみ)

bode関数を利用したゲイン線図のスクリプト

num = [1]; %伝達関数の分子
den = [1 0.5 1];   %伝達関数の分母
minp = -2; %周波数下限のベキ数
maxp = 2;  %周波数上限のベキ数
n = 100;   %周波数分割数

omega = logspace(minp, maxp, n); %周波数ベクトルの生成
sys = tf(num, den); %伝達関数からシステムマトリックスへの変換

[gain phase omega] = bode(sys, omega);
gain = squeeze(gain); % 1次元ベクトルへの変換 gain=gain(:); でもよい
gaindB = 20*log10(gain);

semilogx(omega, gaindB);
xlabel('Frequency (rad/s)');
ylabel('Gain (dB)');
grid( 'on');

周波数応答(ベクトル軌跡)のプロット

ベクトル軌跡は周波数応答関数を$\omega = 0 \sim \infty$の範囲で複素平面上にそのままプロットしたものである.MATLAB(Octave)では複素数ベクトルをplot関数でプロットすると複素平面上で表されるので,それをそのまま利用できる.次のリストがスクリプトの例であり,実行結果は図vector_loci_2_6である.なお,MATLAB(Octave)にはnyquist関数が用意されており,それを利用してもベクトル軌跡を描くことができる.

\frac{1}{s^2 + 0.5s+ 1}のベクトル軌跡

二次遅れ系のベクトル軌跡を描くスクリプト

num = [1]; %伝達関数の分子
den = [1 0.5 1];   %伝達関数の分母
minp = -2; %周波数下限のベキ数
maxp = 2;  %周波数上限のベキ数
n = 100;   %周波数分割数

omega = logspace(minp, maxp, n);
z = polyval(num, complex(0, omega))./polyval(den, complex(0, omega));

plot(z)
axis([-1.5 1.5 -1.5 1.5]);
axis equal; %グラフの縦横比を1:1に指定
grid('on');

応用編

むだ時間要素のボード線図

むだ時間要素$G(s)=e^{-Ls}$は無理式であるので,分子分母にわけて係数ベクトルで表現することはできない.MATLABのControl SystemToolboxではシステムのむだ時間を指定する方法があるが,Octaveではまだ実装されていないようである.ここでは,むだ時間のボード線図を原理に基づいてプロットすることを考える.

むだ時間要素の周波数応答関数は$G(j\omega)=e^{-Lj\omega}$であり,その絶対値と偏角は以下となる.

$$ \begin{array}{l l l} g &=& |G(j \omega)| = 1 \\ \phi &=& \angle G(j \omega) = -L\omega \end{array} $$

よって,むだ時間要素のボード線図を描くスクリプト例は以下のリストとなる.実行結果も示す. ここで,関数onesは指定したサイズで要素がすべて1の行列を生成する関数であり,関数sizeで周波数ベクトルの大きさを取り出して利用している.

むだ時間要素のボード線図(L=1)

むだ時間要素のボード線図を描くスクリプト

L = 1; %むだ時間 [s]
minp = -2; %周波数下限のベキ数
maxp = 2;  %周波数上限のベキ数
n = 100;   %周波数分割数

omega = logspace(minp, maxp, n); %周波数ベクトルの生成
gain = 1*ones(size(omega)); %omegaと同サイズのベクトルを生成
gaindB = 20*log10(gain);
phase = -L*omega/3.14159*180; 

subplot(2,1,1);
semilogx(omega, gaindB);
axis([0.01, 100, -80, 20]);
xlabel('Frequency (rad/s)');
ylabel('Gain (dB)');
grid('on');

subplot(2,1,2);
semilogx(omega, phase);
axis([0.01, 100, -300, 0]);
xlabel('Frequency (rad/s)');
ylabel('Phase (deg)');
grid('on');

実験結果との重ね合わせの例

実験科目や卒業研究などで,実験装置での実測結果と理論式から求めた理論曲線を比較する場面がよく出てくる.次のリストは実験で求めた周波数応答と理論曲線を重ねて描くスクリプトの例である.

実験による測定結果の周波数,ゲイン,位相をそれぞれomega_e, gaindB_e,phase_eという変数に代入しておき,理論曲線と重ねている.この際,実験データの方はマーカーで示すべきであるので,semilog関数内で線類を指定している.また,関数legendでグラフに凡例を付け加えることができる.実行結果を図に示す.

実験結果と理論曲線を重ねる例

実験結果と理論曲線を重ねて描くスクリプト例

omega_e = [0.01 0.03 0.05 0.1 0.3 0.5 1 3 5 10 30 50 100];
gaindB_e = [0 -2.1 3.1 0.1 0.4 1.0 -0.5 -20.1 -27.5 -45.5 -55.2 -69.1 -78.2];
phase_e = [-0.5 -2.1 -3.1 -5.5 -22.1 -32.1 -90.1 -150 -173 -185 -192 -205 -210];

num = [1]; %伝達関数の分子
den = [1 1 1]; %伝達関数の分母
minp = -2; %周波数下限のベキ数
maxp = 2;  %周波数上限のベキ数
n = 100;   %周波数分割数

omega = logspace(minp, maxp, n);
sys = tf(num, den);

[gain phase omega] = bode(sys, omega);
gain = squeeze(gain);
gaindB = 20*log10(gain);
phase = squeeze(phase);

subplot(2,1,1);
semilogx(omega, gaindB, '-', omega_e, gaindB_e, '*');
xlabel('Frequency (rad/s)');
ylabel('Gain (dB)');
grid( 'on');
legend('Theoretical', 'Experimental');

subplot(2,1,2);
semilogx(omega, phase, '-', omega_e, phase_e, '*');
xlabel('Frequency (rad/s)');
ylabel('Phase (deg)');
grid('on');
legend('Theoretical', 'Experimental');