編集の要約なし
編集の要約なし |
編集の要約なし |
||
(2人の利用者による、間の15版が非表示) | |||
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> | ||
38行目: | 38行目: | ||
オイラー法による近似をグラフで示すと図のようになります. | オイラー法による近似をグラフで示すと図のようになります. | ||
[[ファイル:Euler.png]] | |||
次のサンプリング時刻における<math>x(t)</math>の値,すなわち<math>x_{k+1}</math>を接線で一次近似して求めていると解釈できます.もとの微分方程式から,接線の傾きが<math>f(x_k, t_k)</math>ですから,それに<math>\Delta t</math>を乗じた物を<math>x_{k}</math>に加えることで<math>x_{k+1}</math>を求めています.図をみると誤差(error)が大きくて,本当に使えるのか心配ですが,<math>\Delta t</math>をできるだけ小さくすれば誤差を少なくできます.しかし,<math>\Delta t</math>を小さくすると,ある期間(例えば10秒間とか)のシミュレーションをするのにそれだけ計算回数が増えることになるので注意が必要です. | 次のサンプリング時刻における<math>x(t)</math>の値,すなわち<math>x_{k+1}</math>を接線で一次近似して求めていると解釈できます.もとの微分方程式から,接線の傾きが<math>f(x_k, t_k)</math>ですから,それに<math>\Delta t</math>を乗じた物を<math>x_{k}</math>に加えることで<math>x_{k+1}</math>を求めています.図をみると誤差(error)が大きくて,本当に使えるのか心配ですが,<math>\Delta t</math>をできるだけ小さくすれば誤差を少なくできます.しかし,<math>\Delta t</math>を小さくすると,ある期間(例えば10秒間とか)のシミュレーションをするのにそれだけ計算回数が増えることになるので注意が必要です. | ||
44行目: | 44行目: | ||
図では<math>x_k</math>が真値と一致している場合で書きましたが,実際はそうとも限りません.初期値<math>x_0</math>からスタートして,オイラー法を繰り返し使った近似計算は次図のようなイメージになります.<math>\Delta t</math>をできるだけ小さくとって,真値から離れすぎないようにする必要があります.なお,この図では近似のための接線を真値のグラフ<math>x(t)</math>に接するように引いていますが,<math>x_k</math>が真値からずれている場合は図の通りにはなりません.あくまでイメージです. | 図では<math>x_k</math>が真値と一致している場合で書きましたが,実際はそうとも限りません.初期値<math>x_0</math>からスタートして,オイラー法を繰り返し使った近似計算は次図のようなイメージになります.<math>\Delta t</math>をできるだけ小さくとって,真値から離れすぎないようにする必要があります.なお,この図では近似のための接線を真値のグラフ<math>x(t)</math>に接するように引いていますが,<math>x_k</math>が真値からずれている場合は図の通りにはなりません.あくまでイメージです. | ||
[[ファイル:Euler2.png]] | |||
=== Octaveによる実現 === | === Octaveによる実現 === | ||
55行目: | 54行目: | ||
オイラー法を適用すれば | オイラー法を適用すれば | ||
<math>x_{k+1} = x_k + (-x_k + 1) \, \Delta t</math> | |||
とかけるので,繰り返しこの計算を行うOctaveスクリプトの例は以下のようになります. | とかけるので,繰り返しこの計算を行うOctaveスクリプトの例は以下のようになります. | ||
73行目: | 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'); | ||
81行目: | 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>として繰り上げておくことが必要です. | ||
== 高次システムのシミュレーション == | == 高次システムのシミュレーション == | ||
93行目: | 91行目: | ||
<math>\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</math> | <math>\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</math> | ||
[[ファイル:pendulum.png]] | |||
このような2次のシステムの場合,次のような2個の変数を新たに定義します.これらを状態変数と呼びます. | このような2次のシステムの場合,次のような2個の変数を新たに定義します.これらを状態変数と呼びます. | ||
111行目: | 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> | ||
143行目: | 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)'); | ||
== ルンゲ・クッタ法によるシミュレーション == | == ルンゲ・クッタ法によるシミュレーション == | ||
177行目: | 174行目: | ||
[[ファイル: Rungehoine.png]] | |||
*2次のルンゲ・クッタ法(ホイン法)の原理 | *2次のルンゲ・クッタ法(ホイン法)の原理 | ||
201行目: | 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 '); | ||
230行目: | 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)'); | ||
243行目: | 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> | |||
277行目: | 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)'); | ||
316行目: | 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)'); | ||
== 線形システムのシミュレーション == | == 線形システムのシミュレーション == | ||
344行目: | 340行目: | ||
*慣性・粘性・弾性系 | *慣性・粘性・弾性系 | ||
[[ファイル: MDKsystem.png]] | |||
物体の運動方程式は, | 物体の運動方程式は, | ||
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'); |