編集の要約なし
編集の要約なし |
編集の要約なし |
||
| (2人の利用者による、間の14版が非表示) | |||
| 32行目: | 32行目: | ||
<math>\begin{array}{l l l} 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{array}</math> | <math>\begin{array}{l l l} 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{array}</math> | ||
ここで，表記上の簡便性から，<math>x(t_k)</math>を<math>x_k</math>と表記することにし，<math>t_k + \Delta t</math>を<math>t_{k+1}</math>, <math>x(t_k + \Delta t)</math>を | ここで，表記上の簡便性から，<math>x(t_k)</math>を<math>x_k</math>と表記することにし，<math>t_k + \Delta t</math>を<math>t_{k+1}</math>, <math>x(t_k + \Delta t)</math>を<math>x_{k+1}</math>などと表記することにします．この時オイラー法は以下のように表すことができます． | ||
<math>x_{k+1} = x_k + f(x_k, t_k) \, \Delta t</math> | <math>x_{k+1} = x_k + f(x_k, t_k) \, \Delta t</math> | ||
| 54行目: | 54行目: | ||
オイラー法を適用すれば | オイラー法を適用すれば | ||
<math>x_{k+1} = x_k + (-x_k + 1) \, \Delta t</math> | |||
とかけるので，繰り返しこの計算を行うOctaveスクリプトの例は以下のようになります． | とかけるので，繰り返しこの計算を行うOctaveスクリプトの例は以下のようになります． | ||
| 72行目: | 72行目: | ||
endfor | endfor | ||
plot(tim, x); | |||
grid on; | grid on; | ||
axis([0 10 0 1.5]); | axis([0 10 0 1.5]); | ||
xlabel('time (s)'); | xlabel('time (s)'); | ||
ylabel('x'); | ylabel('x'); | ||
| 80行目: | 80行目: | ||
Octaveにおける変数名の都合上，<math>x_k</math>をx_k，<math>x_{k+1}</math>をx_k1としています． | Octaveにおける変数名の都合上，<math>x_k</math>をx_k，<math>x_{k+1}</math>をx_k1としています． | ||
また，シミュレーション結果と時刻を格納する変数をそれぞれx，timとして用意し，<math>x_k</math>と<math>t</math>を行の末尾に追加しています．また，forループで繰り返す際に，次回の計算のために今回の<math>x_{k+1}</math>を次回の<math>x_k</math>として繰り上げておくことが必要です． | また，シミュレーション結果と時刻を格納する変数をそれぞれx，timとして用意し，<math>x_k</math>と<math>t</math>を行の末尾に追加しています．また，forループで繰り返す際に，次回の計算のために今回の<math>x_{k+1}</math>を次回の<math>x_k</math>として繰り上げておくことが必要です． | ||
== 高次システムのシミュレーション == | == 高次システムのシミュレーション == | ||
| 110行目: | 109行目: | ||
<math>\frac{d \boldsymbol{x}(t)}{dt} = \boldsymbol{f}(\boldsymbol{x}(t), t) \quad , \ \boldsymbol{x}(0) = \boldsymbol{x}_0</math> | <math>\frac{d \boldsymbol{x}(t)}{dt} = \boldsymbol{f}(\boldsymbol{x}(t), t) \quad , \ \boldsymbol{x}(0) = \boldsymbol{x}_0</math> | ||
ここで， | ここで，<math>\boldsymbol{f}(\boldsymbol{x}(t), t)</math>は以下となります． | ||
<math>\boldsymbol{f}(\boldsymbol{x}(t), t) = \left(\begin{array}{c}x_2(t)\\- \frac{g}{l} \sin x_1(t) - \frac{c}{m} x_2(t) \end{array}\right)</math> | <math>\boldsymbol{f}(\boldsymbol{x}(t), t) = \left(\begin{array}{c}x_2(t)\\- \frac{g}{l} \sin x_1(t) - \frac{c}{m} x_2(t) \end{array}\right)</math> | ||
| 142行目: | 141行目: | ||
endfor | endfor | ||
plot(tim, x(1, :)); % x_1だけをグラフ表示 | |||
grid on; | grid on; | ||
axis([0 10 -1 1]); | axis([0 10 -1 1]); | ||
xlabel('time (s)'); | xlabel('time (s)'); | ||
ylabel('theta (rad)'); | ylabel('theta (rad)'); | ||
== ルンゲ・クッタ法によるシミュレーション == | == ルンゲ・クッタ法によるシミュレーション == | ||
| 200行目: | 198行目: | ||
endfor | endfor | ||
plot(tim, x); | |||
grid on; | grid on; | ||
axis([0 10 0 1.5]); | axis([0 10 0 1.5]); | ||
xlabel('time (s)'); | xlabel('time (s)'); | ||
ylabel('x '); | ylabel('x '); | ||
| 229行目: | 227行目: | ||
endfor | endfor | ||
plot(tim, x(1, :)); % x_1だけをグラフ表示 | |||
grid on; | grid on; | ||
axis([0 10 -1 1]); | axis([0 10 -1 1]); | ||
xlabel('time (s)'); | xlabel('time (s)'); | ||
ylabel('theta (rad)'); | ylabel('theta (rad)'); | ||
| 242行目: | 240行目: | ||
4次のルンゲ・クッタ法の計算手順は以下のようになります． | 4次のルンゲ・クッタ法の計算手順は以下のようになります． | ||
<math>\begin{array}{l l l}\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{array}</math> | |||
| 276行目: | 274行目: | ||
endfor | endfor | ||
plot(tim, x(1, :)); % x_1だけをグラフ表示 | |||
grid on; | grid on; | ||
axis([0 10 -1 1]); | axis([0 10 -1 1]); | ||
xlabel('time (s)'); | xlabel('time (s)'); | ||
ylabel('theta (rad)'); | ylabel('theta (rad)'); | ||
| 315行目: | 313行目: | ||
endfor | endfor | ||
plot(tim, x(1, :)); % x_1だけをグラフ表示 | |||
grid on; | grid on; | ||
axis([0 10 -1 1]); | axis([0 10 -1 1]); | ||
xlabel('time (s)'); | xlabel('time (s)'); | ||
ylabel('theta (rad)'); | ylabel('theta (rad)'); | ||
== 線形システムのシミュレーション == | == 線形システムのシミュレーション == | ||
| 362行目: | 359行目: | ||
<math>\frac{d}{dt} \left(\begin{array}{c}x_1(t) \\x_2(t)\end{array}\right) =\left(\begin{array}{cc}0 & 1 \\-\frac{K}{M} & -\frac{D}{M}\end{array}\right)\left(\begin{array}{c}x_1(t) \\x_2(t)\end{array}\right)+ \left(\begin{array}{c}0 \\\frac{1}{M}\end{array}\right)u(t)</math> | <math>\frac{d}{dt} \left(\begin{array}{c}x_1(t) \\x_2(t)\end{array}\right) =\left(\begin{array}{cc}0 & 1 \\-\frac{K}{M} & -\frac{D}{M}\end{array}\right)\left(\begin{array}{c}x_1(t) \\x_2(t)\end{array}\right)+ \left(\begin{array}{c}0 \\\frac{1}{M}\end{array}\right)u(t)</math> | ||
状態変数を組み合わせたものを状態ベクトル | 状態変数を組み合わせたものを状態ベクトル<math>\boldsymbol{x}(t)</math>とすれば， | ||
<math>\frac{d\boldsymbol{x}(t)}{dt} =\left(\begin{array}{cc}0 & 1 \\-\frac{K}{M} & -\frac{D}{M}\end{array}\right)\boldsymbol{x}(t)+ \left(\begin{array}{c}0 \\\frac{1}{M}\end{array}\right)u(t)</math> | <math>\frac{d\boldsymbol{x}(t)}{dt} =\left(\begin{array}{cc}0 & 1 \\-\frac{K}{M} & -\frac{D}{M}\end{array}\right)\boldsymbol{x}(t)+ \left(\begin{array}{c}0 \\\frac{1}{M}\end{array}\right)u(t)</math> | ||
となり，これが状態方程式になります．出力は，状態量 | となり，これが状態方程式になります．出力は，状態量<math>x_1(t)</math>そのものですから， | ||
<math>y(t) = \left(\begin{array}{cc}1 & 0\end{array}\right)\boldsymbol{x}(t)</math> | <math>y(t) = \left(\begin{array}{cc}1 & 0\end{array}\right)\boldsymbol{x}(t)</math> | ||
| 375行目: | 372行目: | ||
<math>\boldsymbol{A} = \left(\begin{array}{cc}0 & 1 \\-\frac{K}{M} & -\frac{D}{M}\end{array}\right), \,\boldsymbol{B} = \left(\begin{array}{c}0 \\\frac{1}{M}\end{array}\right), \, \boldsymbol{C} = \left(\begin{array}{cc}1 & 0\end{array}\right), \, \boldsymbol{D} = 0</math> | <math>\boldsymbol{A} = \left(\begin{array}{cc}0 & 1 \\-\frac{K}{M} & -\frac{D}{M}\end{array}\right), \,\boldsymbol{B} = \left(\begin{array}{c}0 \\\frac{1}{M}\end{array}\right), \, \boldsymbol{C} = \left(\begin{array}{cc}1 & 0\end{array}\right), \, \boldsymbol{D} = 0</math> | ||
=== オイラー法によるシミュレーション === | === オイラー法によるシミュレーション === | ||
| 418行目: | 414行目: | ||
endfor | endfor | ||
plot(tim, y); % yをグラフ表示 | |||
grid on; | grid on; | ||
axis([0 10 -1 1]); | axis([0 10 -1 1]); | ||
xlabel('time (s)'); | xlabel('time (s)'); | ||
ylabel('y'); | ylabel('y'); | ||
| 483行目: | 479行目: | ||
endfor | endfor | ||
plot(tim, y); % yをグラフ表示 | |||
grid on; | grid on; | ||
axis([0 10 -1 1]); | axis([0 10 -1 1]); | ||
xlabel('time (s)'); | xlabel('time (s)'); | ||
ylabel('y '); | ylabel('y '); | ||
| 500行目: | 496行目: | ||
以下のような常微分方程式の数値解をlsode関数で求めるための基本的な手順は次のようになります． | 以下のような常微分方程式の数値解をlsode関数で求めるための基本的な手順は次のようになります． | ||
<math>\frac{d \boldsymbol{x}(t)}{dt} = \boldsymbol{f}(\boldsymbol{x}(t), t) \quad , \ \boldsymbol{x}(0) = \boldsymbol{x}_0</math> | |||
*1. まず，微分方程式を適当な関数名（例えばfnc）でOctaveのユーザ関数として定義します．もし解きたい微分方程式が2階以上なら，関数の引数xと戻り値xdotはそれぞれベクトルになります． | *1. まず，微分方程式を適当な関数名（例えばfnc）でOctaveのユーザ関数として定義します．もし解きたい微分方程式が2階以上なら，関数の引数xと戻り値xdotはそれぞれベクトルになります． | ||
| 539行目: | 535行目: | ||
x = lsode("fnc", x_0, t); | x = lsode("fnc", x_0, t); | ||
plot(t, x(:, 1)); % thetaすなわちx(1)をグラフ表示．xは列方向が時間軸になっているの注意 | |||
grid on; | grid on; | ||
axis([0 10 -1 1]); | axis([0 10 -1 1]); | ||
xlabel('time (s)'); | xlabel('time (s)'); | ||
ylabel('theta '); | ylabel('theta '); | ||
なお，このスクリプト例では関数定義の中で質量 | なお，このスクリプト例では関数定義の中で質量<math>m</math>などの定数を定義しています．これらの定数は関数の外側で定義できた方が便利な場合が多いです．通常，関数内部で宣言した変数はローカル変数となり，外部とは名前が同じでも内容は異なるということになってしまうので，グローバルな変数として宣言する必要があります．それにはglobal命令を使います．その際に，関数内部と外部の両方で宣言する必要があるので注意が必要です． | ||
*lsode関数を使った単振り子のシミュレーション2 | *lsode関数を使った単振り子のシミュレーション2 | ||
| 567行目: | 563行目: | ||
x = lsode("fnc", x_0, t); | x = lsode("fnc", x_0, t); | ||
plot(t, x(:, 1)); % thetaすなわちx(1)をグラフ表示．xは列方向が時間軸になっているの注意 | |||
grid on; | grid on; | ||
axis([0 10 -1 1]); | axis([0 10 -1 1]); | ||
xlabel('time (s)'); | xlabel('time (s)'); | ||
ylabel('theta '); | ylabel('theta '); | ||
=== lsode関数による線形システムのシミュレーション === | === lsode関数による線形システムのシミュレーション === | ||
| 580行目: | 575行目: | ||
<math>\begin{array}{l l l}\frac{d \boldsymbol{x}(t)}{dt} &=& \boldsymbol{A} \boldsymbol{x}(t) + \boldsymbol{B} \boldsymbol{u}(t) \\\boldsymbol{y}(t) &=& \boldsymbol{C} \boldsymbol{x}(t) + \boldsymbol{D} \boldsymbol{u}(t)\end{array}</math> | <math>\begin{array}{l l l}\frac{d \boldsymbol{x}(t)}{dt} &=& \boldsymbol{A} \boldsymbol{x}(t) + \boldsymbol{B} \boldsymbol{u}(t) \\\boldsymbol{y}(t) &=& \boldsymbol{C} \boldsymbol{x}(t) + \boldsymbol{D} \boldsymbol{u}(t)\end{array}</math> | ||
lsode関数で計算するのは状態方程式のほうだけです．その結果得られた状態変数ベクトル<math>\boldsymbol{x}</math>を用いて出力方程式に基づいて出力ベクトル | lsode関数で計算するのは状態方程式のほうだけです．その結果得られた状態変数ベクトル<math>\boldsymbol{x}</math>を用いて出力方程式に基づいて出力ベクトル<math>\boldsymbol{y}</math>を計算します． | ||
ここでは前述の慣性・粘性・弾性系の強制振動をlsode関数を用いてシミュレーションするスクリプト例を示します．このスクリプトで注意すべきなのは，入力<math>u(t)</math>の取り扱いです．<math>u(t)</math>はlsodeの内部では変数tの代数式として記述しますが，lsodeで得られた状態ベクトルを用いて出力yを計算する際には時間軸ベクトルtに対応したベクトルuとして必要になります．Octaveではsin関数などが引数にベクトルを受けると，要素ごとに関数をかけて同じサイズのベクトルを返しますのでそれを利用しています．<math>u(t)</math>の内容によっては，要素ごとの積算 .* や要素ごとの除算 ./ といった演算子が必要になるので注意です． | ここでは前述の慣性・粘性・弾性系の強制振動をlsode関数を用いてシミュレーションするスクリプト例を示します．このスクリプトで注意すべきなのは，入力<math>u(t)</math>の取り扱いです．<math>u(t)</math>はlsodeの内部では変数tの代数式として記述しますが，lsodeで得られた状態ベクトルを用いて出力yを計算する際には時間軸ベクトルtに対応したベクトルuとして必要になります．Octaveではsin関数などが引数にベクトルを受けると，要素ごとに関数をかけて同じサイズのベクトルを返しますのでそれを利用しています．<math>u(t)</math>の内容によっては，要素ごとの積算 .* や要素ごとの除算 ./ といった演算子が必要になるので注意です． | ||
| 613行目: | 608行目: | ||
y = C * x '+ D * u; %出力方程式に基づいてyの計算 | y = C * x '+ D * u; %出力方程式に基づいてyの計算 | ||
plot(t, y); % yをグラフ表示 | |||
grid on; | grid on; | ||
axis([0 10 -1 1]); | axis([0 10 -1 1]); | ||
xlabel('time (s)'); | xlabel('time (s)'); | ||
ylabel('y '); | ylabel('y '); | ||
== 線形離散時間システムとしてのシミュレーション == | |||
これまで常微分方程式の数値解法に基づいて動的システムのシミュレーション方法を考えてきましたが，シミュレーション対象が線形システムで，入力がサンプリング期間中一定値と考えてよい場合，線形離散時間システムとして扱うことで簡単にシミュレーションを行うことができます． | これまで常微分方程式の数値解法に基づいて動的システムのシミュレーション方法を考えてきましたが，シミュレーション対象が線形システムで，入力がサンプリング期間中一定値と考えてよい場合，線形離散時間システムとして扱うことで簡単にシミュレーションを行うことができます． | ||
| 634行目: | 628行目: | ||
です． | です． | ||
ここで，<math>\boldsymbol{\Phi} | ここで，<math>\boldsymbol{\Phi}</math>と<math>\boldsymbol{\Gamma}</math>が求まれば，離散時間システムの式を直接演算するだけで行えるので簡単です． | ||
連続時間システムから離散時間システムへの変換には，Octaveに用意されているc2d関数が利用できます．この関数はサンプル・ホールド付きのz変換か双一次変換のいずれかで離散化することができます（デフォルトはサンプル・ホールド）．シミュレーションにはサンプル・ホールド付きのz変換を使います． | 連続時間システムから離散時間システムへの変換には，Octaveに用意されているc2d関数が利用できます．この関数はサンプル・ホールド付きのz変換か双一次変換のいずれかで離散化することができます（デフォルトはサンプル・ホールド）．シミュレーションにはサンプル・ホールド付きのz変換を使います． | ||
| 640行目: | 634行目: | ||
次のスクリプトは，先の慣性・粘性・弾性系の強制振動を離散時間システムに変換した上でシミュレーションするものです．まず，ss関数で4つの行列による状態空間表現をOctave上のシステム行列表現に変換し，その上でc2d関数で離散化しています．シミュレーションの計算自体は離散時間システムの式を繰り返し演算しているだけです． | 次のスクリプトは，先の慣性・粘性・弾性系の強制振動を離散時間システムに変換した上でシミュレーションするものです．まず，ss関数で4つの行列による状態空間表現をOctave上のシステム行列表現に変換し，その上でc2d関数で離散化しています．シミュレーションの計算自体は離散時間システムの式を繰り返し演算しているだけです． | ||
=== 離散時間システムに変換した慣性・粘性・弾性系のシミュレーション === | |||
dt = 0.01; % サンプリングタイム | dt = 0.01; % サンプリングタイム | ||
| 675行目: | 669行目: | ||
endfor | endfor | ||
plot(tim, y); % yをグラフ表示 | |||
grid on; | grid on; | ||
axis([0 10 -1 1]); | axis([0 10 -1 1]); | ||
xlabel('time (s)'); | xlabel('time (s)'); | ||
ylabel('y'); | ylabel('y'); | ||