編集の要約なし
編集の要約なし |
編集の要約なし |
||
| 79行目: | 79行目: | ||
ylabel('x'); | ylabel('x'); | ||
Octaveにおける変数名の都合上， | Octaveにおける変数名の都合上，<math>x_k</math>をx_k，<math>x_{k+1</math>をx_k1としています． | ||
また，シミュレーション結果と時刻を格納する変数をそれぞれx，timとして用意し， | また，シミュレーション結果と時刻を格納する変数をそれぞれx，timとして用意し，<math>x_k</math>と<math>t</math>を行の末尾に追加しています．また，forループで繰り返す際に，次回の計算のために今回の<math>x_{k+1}</math>を次回の<math>x_k</math>として繰り上げておくことが必要です． | ||
| 89行目: | 89行目: | ||
これまでは1階(次)の微分方程式を扱ってきましたが，現実には2次以上のシステムのシミュレーションが必要になることが多いと思います．その場合は，変数変換をすることで連立多元1階常微分方程式に変形して考えます． | これまでは1階(次)の微分方程式を扱ってきましたが，現実には2次以上のシステムのシミュレーションが必要になることが多いと思います．その場合は，変数変換をすることで連立多元1階常微分方程式に変形して考えます． | ||
例えば，次のような微分方程式を考えます．これは，長さ | 例えば，次のような微分方程式を考えます．これは，長さ<math>l</math>の軽い棒の先に質量<math>m</math>の重りがついた単振り子を初期角度<math>\theta_0</math>の位置から初速度0で放す状況を表しています（下図参照）．ここで振り子は空気その他から速度に比例した抵抗を受けるものとし，その係数を<math>c</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> | |||
#ref(pendulum.png) | #ref(pendulum.png) | ||
| 97行目: | 97行目: | ||
このような2次のシステムの場合，次のような2個の変数を新たに定義します．これらを状態変数と呼びます． | このような2次のシステムの場合，次のような2個の変数を新たに定義します．これらを状態変数と呼びます． | ||
<math>\begin{array}{l l l}x_1(t) &=& \theta(t)\\x_2(t) &=& \frac{d \theta(t)}{dt}\end{array}</math> | |||
状態変数を用いてもとの微分方程式を表すと以下のように1階の微分方程式が2本連立したものになります． | 状態変数を用いてもとの微分方程式を表すと以下のように1階の微分方程式が2本連立したものになります． | ||
<math>\begin{array}{l l l}\frac{d x_1(t)}{dt} &=& x_2(t) \quad , x_1(0) = \theta_0\\\frac{d x_2(t)}{dt} &=& - \frac{g}{l} \sin x_1(t) - \frac{c}{m} x_2(t) \quad , x_2(0) = 0\end{array}</math> | |||
さらに，2つの状態変数を組にして， | さらに，2つの状態変数を組にして， | ||
<math>\boldsymbol{x}(t) = \left(\begin{array}{c}x_1(t)\\x_2(t)\end{array}\right) \quad , \ \boldsymbol{x}_0 = \left(\begin{array}{c}\theta_0\\0\end{array}\right)</math> | |||
のようにベクトル | のようにベクトル<math>\boldsymbol{x}(t)</math>で表現することで，最初の1階常微分方程式と同じ見た目になります． | ||
<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(\boldsymbol{f}(\boldsymbol{x}(t), t));は以下となります． | ||
<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{x}_{k+1} = \boldsymbol{x}_k + \boldsymbol{f}(\boldsymbol{x}_k, t_k) \, \Delta t</math> | |||
この状態変数表現に基づいたOctaveによるシミュレーション・スクリプトは次のようになります． | この状態変数表現に基づいたOctaveによるシミュレーション・スクリプトは次のようになります． | ||
基本的には1次のシステムに対するオイラー法のスクリプトと同じですが，x_k やx_k1が2行×1列の列ベクトルになっていることが違いです．それに伴い，結果を格納する変数xが最終的に2行×サンプリング点数の行列になります．よって，シミュレーション結果として | 基本的には1次のシステムに対するオイラー法のスクリプトと同じですが，x_k やx_k1が2行×1列の列ベクトルになっていることが違いです．それに伴い，結果を格納する変数xが最終的に2行×サンプリング点数の行列になります．よって，シミュレーション結果として<math>\theta</math>だけをグラフ表示するには1行目だけを取り出す必要があり，そのため，plot関数の引数をx(1, : )としています． | ||
なお，この例ではサンプリングタイム | なお，この例ではサンプリングタイム<math>\Delta t</math>を0.01秒としていますが，これを0.05秒にするとシミュレーション結果が大きく異なることに気がつくと思います．さらに<math>\Delta t</math>を大きくしていくと，それは誤差が大きくなるというようなレベルではなく，全く異なったシミュレーション結果になってしまいます．このようにサンプリングタイムの選択は重要で，基本的には十分小さくとる必要があります． | ||
*オイラー法による単振り子のシミュレーション | *オイラー法による単振り子のシミュレーション | ||
| 154行目: | 154行目: | ||
=== 2次のルンゲ・クッタ法（ホイン法） === | === 2次のルンゲ・クッタ法（ホイン法） === | ||
オイラー法は | オイラー法は<math>t_k</math>における導関数の値<math>f(x_k, t_k)</math>を使って<math>x_{k+1}</math>を近似していました．ここでは．より精度を高めるために，複数点における導関数の値を利用することを考えます． | ||
先の図のように | 先の図のように<math>t_k</math>における接線の傾き<math>f(x_k, t_k)</math>だけを考えるのではなく，<math>t_{k+1}</math>における傾き<math>f(x_{k+1}, t_{k+1})</math>も利用して，それらの平均をその区間での傾きとして利用すれば，精度が高まりそうです．その場合，<math>x_{k+1}</math>は以下のように考えられます． | ||
<math>x_{k+1} = x_k + \frac{1}{2} ( f(x_k, t_k) + f(x_{k+1}, t_{k+1}) ) \Delta t</math> | |||
しかし， | しかし，<math>x_{k+1}</math>を求める計算に，<math>x_{k+1}</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 + \frac{1}{2} ( f(x_k, t_k) + f(x_k + f(x_k, t_k) \Delta t, t_{k+1}) ) \Delta t</math> | |||
とかけます．さらに，逐次計算しやすいように | とかけます．さらに，逐次計算しやすいように<math>k_1</math>, <math>k_2</math>という変数を定義して整理すれば， | ||
<math>\begin{array}{l l l}k_1 &=& f(x_k, t_k) \Delta t \\k_2 &=& f(x_k + k_1, t_{k+1}) \Delta t \\x_{k+1} &=& x_k + \frac{1}{2} ( k_1 + k_2 )\end{array}</math> | |||
となります．これは2次のルンゲ・クッタ法の一種で，ホイン法または修正オイラー法と呼ばれるものです．2次のルンゲ・クッタ法としては，導関数を取る点として区間の中点を考えた中点法と呼ばれるものもあります． | となります．これは2次のルンゲ・クッタ法の一種で，ホイン法または修正オイラー法と呼ばれるものです．2次のルンゲ・クッタ法としては，導関数を取る点として区間の中点を考えた中点法と呼ばれるものもあります． | ||
| 181行目: | 181行目: | ||
*2次のルンゲ・クッタ法（ホイン法）の原理 | *2次のルンゲ・クッタ法（ホイン法）の原理 | ||
前節でオイラー法を用いてシミュレーションした1階の微分方程式と，2階の微分方程式をそれぞれホイン法でシミュレーションするOctaveスクリプトは以下のようになります．2階の微分方程式の場合は | 前節でオイラー法を用いてシミュレーションした1階の微分方程式と，2階の微分方程式をそれぞれホイン法でシミュレーションするOctaveスクリプトは以下のようになります．2階の微分方程式の場合は<math>k_1</math>, <math>k_</math>もベクトルになることに注意してください． | ||
また，オイラー法に比べてサンプリングタイムを大きくとっても精度よくシミュレーションできることを確かめてみてください． | また，オイラー法に比べてサンプリングタイムを大きくとっても精度よくシミュレーションできることを確かめてみてください． | ||