編集の要約なし
編集の要約なし |
編集の要約なし |
||
| 79行目: | 79行目: | ||
ylabel('x'); | ylabel('x'); | ||
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>として繰り上げておくことが必要です． | ||
| 360行目: | 360行目: | ||
これを行列を用いて表現すると以下のようになります． | これを行列を用いて表現すると以下のようになります． | ||
<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(\boldsymbol{x}(t));とすれば， | ||
<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(x_1(t));そのものですから， | ||
<math>y(t) = \left(\begin{array}{cc}1 & 0\end{array}\right)\boldsymbol{x}(t)</math> | |||
となります． | となります． | ||
以上から，状態空間表現の行列 | 以上から，状態空間表現の行列<math>\boldsymbol{A},\boldsymbol{B},\boldsymbol{C},\boldsymbol{D}</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> | |||
| 383行目: | 383行目: | ||
状態方程式に対するオイラー法の演算は以下のように書き表せます． | 状態方程式に対するオイラー法の演算は以下のように書き表せます． | ||
<math>\begin{array}{l l l}\boldsymbol{x}_{k+1} &=& \boldsymbol{x}_k + ( \boldsymbol{A} \boldsymbol{x}_k + \boldsymbol{B} \boldsymbol{u}_k ) \Delta t \\\boldsymbol{y}_{k} &=& \boldsymbol{C} \boldsymbol{x}_k + \boldsymbol{D} \boldsymbol{u}_k \end{array}</math> | |||
この式に基づいて先の慣性・粘性・弾性系のシミュレーションを行うOctaveスクリプトの例を示します．入力 | この式に基づいて先の慣性・粘性・弾性系のシミュレーションを行うOctaveスクリプトの例を示します．入力<math>\boldsymbol{u}(t)</math>を正弦波とした強制振動を想定しています． | ||
*オイラー法による慣性・粘性・弾性系のシミュレーション | *オイラー法による慣性・粘性・弾性系のシミュレーション | ||
| 428行目: | 428行目: | ||
状態空間表現された線形システムを(4次の)ルンゲ・クッタ法でシミュレーションする演算式は以下のようになります． | 状態空間表現された線形システムを(4次の)ルンゲ・クッタ法でシミュレーションする演算式は以下のようになります． | ||
<math>\begin{array}{l l l}\boldsymbol{k}_1 &=& \left( \boldsymbol{A} \boldsymbol{x}_k + \boldsymbol{B} \boldsymbol{u}_k \right) \Delta t \\\boldsymbol{k}_2 &=& \left( \boldsymbol{A} (\boldsymbol{x}_k + \frac{\boldsymbol{k}_1}{2}) + \boldsymbol{B} \boldsymbol{u}(t_k + \frac{\Delta t}{2}) \right) \Delta t \\\boldsymbol{k}_3 &=& \left( \boldsymbol{A} (\boldsymbol{x}_k + \frac{\boldsymbol{k}_2}{2}) + \boldsymbol{B} \boldsymbol{u}(t_k + \frac{\Delta t}{2}) \right) \Delta t \\\boldsymbol{k}_4 &=& \left( \boldsymbol{A} (\boldsymbol{x}_k + \boldsymbol{k}_3) + \boldsymbol{B} \boldsymbol{u}_{k+1} \right) \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)\\\boldsymbol{y}_{k} &=& \boldsymbol{C} \boldsymbol{x}_k + \boldsymbol{D} \boldsymbol{u}_k \end{array}</math> | |||
ルンゲ・クッタ法の場合，時刻 | ルンゲ・クッタ法の場合，時刻<math>t_k</math>だけでなく，<math>t_k + \frac{\Delta t}{2}</math>や<math>t_{k+1}</math>での入力<math>\boldsymbol{u}(t)</math>が必要になることに注意します． | ||
ルンゲ・クッタ法で慣性・粘性・弾性系のシミュレーションを行うOctaveスクリプトの例は以下の通りです． | ルンゲ・クッタ法で慣性・粘性・弾性系のシミュレーションを行うOctaveスクリプトの例は以下の通りです． | ||
| 492行目: | 492行目: | ||
=== 多段解法 === | === 多段解法 === | ||
これまで微分方程式の数値解法としてオイラー法とルンゲ・クッタ法について述べましたが，両手法とも | これまで微分方程式の数値解法としてオイラー法とルンゲ・クッタ法について述べましたが，両手法とも<math>x_{k+1}</math>を求めるために状態量としては<math>x_k</math>だけを利用する方法でした（一段階法）．しかし，<math>x_{k-1}</math>や<math>x_{k-2}</math>といった，より過去の情報も利用する多段階法という手法があります．複数の過去の情報を利用することで計算精度を上げることができるだけでなく，数値解法の安定性の点でも優れることがわかっています．すなわち，オイラー法やルンゲ・クッタ法では，サンプリング周期<math>\Delta t</math>を大きくとりすぎるとシミュレーション結果が発散してしまうことがありましたが，多段階法では発散しにくいのです．この性質は，時間のかかるシミュレーションをする前に<math>\Delta t</math>を大きめにとってラフなシミュレーションをすることができるので便利です． | ||
=== lsode関数の使用方法 === | === lsode関数の使用方法 === | ||