本演習では，走行のシンプルモデルとして知られるSpring-Loaded Inverted Pendulum (SLIP)モデルの動力学シミュレーションを構築し，走行の力学について考えます．シミュレータはPythonを用いて自分で構築してもらいます． 第1週ではまず，動力学シミュレーションを作る練習として，単純なバネマスダンパ系の運動について考えましょう．
まずは講義形式で動力学と数値計算について簡単に復習を行います．その後，各自のPCを用いて演習を行います．
Pythonをあまり使ったことのない方は，まずこちらを確認してください．Python環境の準備についても記載しています．
質点$m$がバネとダンパに繋がれた系を考える．バネとダンパの係数はそれぞれ$k$と$c$とする．バネの平衡点から測定した質点の位置を$x$として，系の運動方程式を求め，その運動を厳密解と数値解を用いて観察しよう．なお，系に外力は働かないものとする．また，$\dot{*}$は変数$$*を時間微分したものを表す．
(1) $x$の2階微分方程式の形で，系の運動方程式を書き下そう．
(2) $v=\dot{x}$として，$\bm{y}=\begin{bmatrix}x\\ v\end{bmatrix}$とするとき，(1)で求めた2階微分方程式を，ベクトル$\bm{y}$に関する1階の微分方程式に書き直そう．
(3) (1)または(2)で導出した微分方程式の解析解$x(t)$を求めよう．初期条件は$x(0)=x_0$，$\dot{x}(0)=0$とする．
(4) Pythonを用いて，オイラー法によって(2)の微分方程式を$t\in[0,3]$ [s]の範囲で数値積分し，その結果を$t$-$x$平面にプロットするプログラムを作ろう．なお，$m=1.0$ [kg]，$k=100$ [N/m]，$c = 10$ [Ns/m]とし，数値計算の時間間隔は$\Delta t=0.05$ [s]とする．Pythonの使い方はこちら．
(5) (3)で求めた厳密解を，(4)で求めた数値解と併せてプロットし，オイラー法の計算精度を確認しよう．
(6) (4)で作成したプログラムをもとにして，4次のルンゲ・クッタ法によって(2)の微分方程式を数値積分するプログラムを作ろう．条件は(4)と同一とする．
(7) (6)で求めた数値解と(3)の厳密解と比較し，オイラー法よりもルンゲ・クッタ法の方が計算精度に優れていることを確認しよう．応用課題：$\Delta t$を横軸に，誤差を縦軸にして，両対数グラフにプロットし，オイラー法とルンゲクッタ法の1ステップでの誤差がそれぞれ$\Delta t^2$と$\Delta t^5$のオーダーになっているかも確認できるとなお良い（$\Delta t$は0.01から0.10くらいまで変化させ，最初の1ステップでの誤差を比べる）．そうでない場合は理由も考察してみよ．
(8) 修正オイラー法は，以下によって数値積分を行うアルゴリズムで，オイラー法よりも精度が優れている．
$$ \begin{align*} \tilde{x}(t+\Delta t) &:=x(t) + f(x(t))\Delta t\\ {x}(t+\Delta t) &\approx x(t)+\frac{1}{2}\left\{f(x(t))+f(\tilde{x}(t+\Delta t))\right\}\Delta t \end{align*} $$
厳密解に対してオイラー法の誤差は$O(\Delta t^2)$であり，修正オイラー法の誤差は$O(\Delta t^3)$であることを証明しよう．（なお，ルンゲ・クッタ法の誤差はさらに小さく$O(\Delta t^5)$である．）
(9) 質点に強制外力$F=A\sin\Omega t$を加えたときの挙動を，解析的にも数値的にも確認しよう．