編集の要約なし
編集の要約なし |
|||
| (2人の利用者による、間の6版が非表示) | |||
| 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'); | ||
| 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)'); | ||
| 198行目: | 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 '); | ||
| 227行目: | 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)'); | ||
| 274行目: | 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)'); | ||
| 313行目: | 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)'); | ||
| 359行目: | 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> | ||
| 372行目: | 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> | ||
=== オイラー法によるシミュレーション === | === オイラー法によるシミュレーション === | ||
| 415行目: | 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'); | ||
| 480行目: | 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 '); | ||
| 536行目: | 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 '); | ||
| 564行目: | 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 '); | ||
| 609行目: | 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 '); | ||
== 線形離散時間システムとしてのシミュレーション == | |||
これまで常微分方程式の数値解法に基づいて動的システムのシミュレーション方法を考えてきましたが，シミュレーション対象が線形システムで，入力がサンプリング期間中一定値と考えてよい場合，線形離散時間システムとして扱うことで簡単にシミュレーションを行うことができます． | これまで常微分方程式の数値解法に基づいて動的システムのシミュレーション方法を考えてきましたが，シミュレーション対象が線形システムで，入力がサンプリング期間中一定値と考えてよい場合，線形離散時間システムとして扱うことで簡単にシミュレーションを行うことができます． | ||
| 635行目: | 634行目: | ||
次のスクリプトは，先の慣性・粘性・弾性系の強制振動を離散時間システムに変換した上でシミュレーションするものです．まず，ss関数で4つの行列による状態空間表現をOctave上のシステム行列表現に変換し，その上でc2d関数で離散化しています．シミュレーションの計算自体は離散時間システムの式を繰り返し演算しているだけです． | 次のスクリプトは，先の慣性・粘性・弾性系の強制振動を離散時間システムに変換した上でシミュレーションするものです．まず，ss関数で4つの行列による状態空間表現をOctave上のシステム行列表現に変換し，その上でc2d関数で離散化しています．シミュレーションの計算自体は離散時間システムの式を繰り返し演算しているだけです． | ||
=== 離散時間システムに変換した慣性・粘性・弾性系のシミュレーション === | |||
dt = 0.01; % サンプリングタイム | dt = 0.01; % サンプリングタイム | ||
| 670行目: | 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'); | ||