編集の要約なし
編集の要約なし |
|||
| (2人の利用者による、間の10版が非表示) | |||
| 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> | ||
| 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'); | ||
| 109行目: | 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> | ||
| 141行目: | 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)'); | ||
== ルンゲ・クッタ法によるシミュレーション == | == ルンゲ・クッタ法によるシミュレーション == | ||
| 199行目: | 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 '); | ||
| 228行目: | 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)'); | ||
| 241行目: | 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> | |||
| 275行目: | 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)'); | ||
| 314行目: | 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)'); | ||
== 線形システムのシミュレーション == | == 線形システムのシミュレーション == | ||
| 361行目: | 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> | ||
| 374行目: | 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> | ||
=== オイラー法によるシミュレーション === | === オイラー法によるシミュレーション === | ||
| 417行目: | 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'); | ||
| 482行目: | 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 '); | ||
| 538行目: | 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 '); | ||
| 566行目: | 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 '); | ||
| 578行目: | 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>の内容によっては，要素ごとの積算 .* や要素ごとの除算 ./ といった演算子が必要になるので注意です． | ||
| 611行目: | 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 '); | ||
== 線形離散時間システムとしてのシミュレーション == | |||
これまで常微分方程式の数値解法に基づいて動的システムのシミュレーション方法を考えてきましたが，シミュレーション対象が線形システムで，入力がサンプリング期間中一定値と考えてよい場合，線形離散時間システムとして扱うことで簡単にシミュレーションを行うことができます． | これまで常微分方程式の数値解法に基づいて動的システムのシミュレーション方法を考えてきましたが，シミュレーション対象が線形システムで，入力がサンプリング期間中一定値と考えてよい場合，線形離散時間システムとして扱うことで簡単にシミュレーションを行うことができます． | ||
| 638行目: | 634行目: | ||
次のスクリプトは，先の慣性・粘性・弾性系の強制振動を離散時間システムに変換した上でシミュレーションするものです．まず，ss関数で4つの行列による状態空間表現をOctave上のシステム行列表現に変換し，その上でc2d関数で離散化しています．シミュレーションの計算自体は離散時間システムの式を繰り返し演算しているだけです． | 次のスクリプトは，先の慣性・粘性・弾性系の強制振動を離散時間システムに変換した上でシミュレーションするものです．まず，ss関数で4つの行列による状態空間表現をOctave上のシステム行列表現に変換し，その上でc2d関数で離散化しています．シミュレーションの計算自体は離散時間システムの式を繰り返し演算しているだけです． | ||
=== 離散時間システムに変換した慣性・粘性・弾性系のシミュレーション === | |||
dt = 0.01; % サンプリングタイム | dt = 0.01; % サンプリングタイム | ||
| 673行目: | 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'); | ||