運動方程式の解法 - 物理学の見つけ方

ニュートンの運動方程式：を解きたい。

式()を数値的に解きたい

前章では、この方程式を導出するところまで行ったので、後は解き方が分かれば、物体の運動を実際に計算できるようになる。

この章では、式()の一般的な解放である数値解法について述べた後、いくつかの例について実際に計算を行う： $運動方程式の数値解法力の例（参考）オイラー法は微分方程式の解に収束する$

3.1運動方程式の数値解法

ほとんどの運動方程式は、解析的には解けないので、数値的に解くことになる。この節では、もっとも単純な解法であるオイラー法()について述べる。（その数学的根拠については3節で示す。）

オイラー法による運動方程式の解法：式()

運動方程式()を数値的に解く方法は、第1章のキャッチボールの場合と同じである。即ち、1次近似によって得られる以下の漸化式：を用いて、ずつ時刻を進めていけばよい。ただし、漸化式の初期値、即ち、での位置と速度は与えられているとする。この式は、微分方程式を解く一般的な方法であるオイラー法（以下の【3.1-注1】）を適用したものである。

【3.1-注1】微分方程式の解法：オイラー法

初期値が与えられた時、1階正規形の微分方程式：（以下の【3.1-注2】の()）を数値的に解いて任意の時刻でのを求めたい。それには、1次近似の式：を漸化式として使って、十分に小さな時間間隔ずつ時刻を進めていけばよい（右図）。この計算手法を、オイラー法という。

補足

オイラー法は単純で分かりやすい解法ではあるが、を相当小さくしなければ誤差が十分減らない。しかし一方、を小さくしすぎると、計算量が増大することに加え、丸め誤差が問題になる。丸め誤差とは、コンピュータ上では実数が有限の精度しか持たないことによる誤差であり、例えば精度が15桁の場合、はに丸められてしまい変化しなくなる（が小さすぎる場合の式()の右辺に対応する）。これらのため実用上は、オイラー法()のようにの1次関数で近似するのではなく、より高次の多項式で近似する方法がとられる。特に、4次のルンゲ・クッタ法と呼ばれるものがよく使われ、本サイトのシミュレーションでも使用しているが、説明は割愛する。

【3.1-注2】1階正規形の微分方程式()

未知の関数に対する以下の形の式を、1階正規形の微分方程式と呼ぶ：（は既知の関数）「1階」はについての1階微分方程式であることを意味し、「正規形」はの形になっていることを意味する。

補足

多くの微分方程式がこの形に変形できる。特に、運動方程式は2階の微分方程式であるが、以下のように1階正規形に変形できる：赤字部分をとおき、最後の式をとおけば、確かに1階正規形()になっている。

力は、時間・位置・速度の関数

オイラー法()で運動方程式を解くには、右辺に現れる力が分かっている必要がある。は、時間の既知の関数でなくとも、未知関数である位置や速度の関数であってもよい：（オイラー法()の右辺ではは既知だからである。）

もう少し詳しく言うと、の依存性では、例えばのように、どのような関数になるか予め分かっている必要がある。の依存性は、例えばバネがそうであったように、その時刻での位置が分かった時点で、そのの値に応じて働く力が決まるというものである。の依存性も同様で、例としては空気抵抗が挙げられる（速度が大きいほど抵抗が大きくなる）。

3.2力の例

この節では、いくつかの例に対して、シミュレーションを行う。単純なものについては解析解も示す。

【例1】一定の力：式()

もっとも単純な例として、が定数の場合、即ち、時間や位置などによらず常に一定の力で押されて（引かれて）いる場合を考える。このを表したのが右図である。右図はベクトル場と呼ばれるもので、「物体をそこに置いた時にどのような力が働くか」が矢印で示されている。こを見れば、物体がどこにあっても上向きの一定の力が働くことが一目瞭然である。

運動方程式は、である。

数値計算を行うと右図のようになる。式()は、前章のキャッチボールの運動方程式とほとんど同じ形なので、解析解が容易に見つかる（以下の【3.2-注1】）。

【3.2-注1】定数の力の場合の解析解：式()

運動方程式()の解析解は、を初期値として、以下のようになる：

【例2】変位に比例する力：式()

定数の次に単純なのは、力がに比例する場合である：（簡単のため軸上の1次元運動とする）。この力をベクトル場で表すと、右図のようになる（上側は、下側はの場合）。

運動方程式はである。

数値計算を行うと右図のようになる。の場合、原点に引きつけようとする力を受けて、振動する（例えばバネや振れ幅の小さい振り子）。逆に、の場合、原点から勢いよく離れていく（例えば逆さまにした振り子の頂点付近）。解析解も知られており、以下の【3.2-注2】のようになる。

式()のような運動方程式は、現実にもよく現れる。それは、釣合いの点の付近である。実際、物体の位置が釣合いの点（＝原点とする）の付近にある場合、力を1次近似すれば、確かにの形になる：はでの微分である：（一方、と書くと時間微分）。空気抵抗など、速度に依存する力は無視できるとする。特に、バネの場合、この近似が比較的広い範囲で成り立つ（以下の【3.2-注3】）。

【3.2-注2】1次の力の場合の解析解：式()

運動方程式()の解析解は、の符号によって異なり、以下のようになる：はそれぞれ初期位置・初期速度であり、は以下で定義される：に現れるは、双曲線関数である：

【3.2-注3】フックの法則

実験によると、「バネの長さの変化」と「バネに加えた力」は、近似的に比例する。即ち、バネに固有の定数（＝バネ定数という）を用いて、以下が成り立つ：この関係式を、フックの法則という（でもよい）。フックの法則はバネ量りに使える。実際、バネ定数をあらかじめ求めておけば、任意の物体を吊り下げたときのを読み取ることにより、その物体に働く重力（＝重さ）が逆算できる。

【例3】惑星に働く力：式()

地球は、太陽の周りを回っている。即ち、太陽から力を受けている。この力は、原点に太陽があるとして、以下の式で表せることが知られている：（は地球の位置、は地球の質量。）は、太陽の方向を向き、大きさは太陽からの距離の乗で減衰する。詳しくは重力理論編の第1章を参照。

数値計算を行うと右図のようになる。

【例4】磁場中の荷電粒子に働く力：式()

電流や磁石の近くを運動する荷電粒子は、速度に依存した力（ローレンツ力）を受けることが知られている。この力は以下のようになる：（右辺は行列とベクトルの積）はに垂直である（との内積をとればゼロになる）。は磁場と呼ばれるベクトルで、位置や時間に依存してもよい。詳しくは電磁力学編の第3章を参照。

一様磁場の場合に数値計算を行うと右図のようになる。

3.3（参考）オイラー法は微分方程式の解に収束する

この節では、オイラー法のステップ幅を小さくしていくと、解に収束することを示す。ここでは解の存在を仮定するが、解が存在すること自体もほぼ同じ条件の下で示せる。

定理：オイラー法は微分方程式の解に収束する

微分方程式が、任意の初期値の下で常に解を持つとする。この時、オイラー法でステップ幅の極限をとったものは、その解に一致する。

ただし、の変化率は有限であるとする。即ち、ある定数（＝傾きの最大値）が存在して、（同一時刻の）任意の異なる2点に対してが成り立つ。これをリプシッツ条件という。リプシッツ条件()のノルムは、のようにとればよい。異なる単位が混ざった2乗和が現れるとそのままでは定義できないが、その場合は、単位を無視して数値部分だけを使えばよい（単位の取り方によって値が変わってしまうが「傾きが有限」という性質には影響しないので問題ない）。

証明

微分方程式()の厳密解を、オイラー法で得られた折れ線関数をとおく（はオイラー法の刻み幅）。示したいのは、の時に両者が一致することである：そのためには、厳密解とオイラー法の誤差の最大値を見積もり（＝以下の式()）、その誤差がの時にゼロになることを示せばよい。

まず、オイラー法の1ステップの間に生じ得る、最大誤差をとおく（右図）。はステップ幅だけの単調増加関数であり、オイラー法の起点には依存しないとする（あらゆる起点を考えた上での最大値を考えている）。は、1次近似で生じる誤差なので、の2次程度の微小量（）：である。

それでは証明に入ろう。時刻からの間に生じ得る最大誤差をとおく。オイラー法の1ステップでの誤差は右図のように評価できるので、とを関係づける以下の漸化式を得る：青色部分は1ステップ前の最大誤差を引き継いだもの、緑色部分は傾きの誤差の影響をリプシッツ条件()（分母をに置き換えたもの）を使って評価したもの、マゼンタ部分は式()である。数式的に式()を導くには、三角不等式を使っての左辺を、式()の各項が出るように分解していけばよい。漸化式()を、初期値の下で解くと $後のために、公式を適用する$ となる（特性方程式を使って解ける）。

式()において、を保ったままの極限を取ればとなる（式()により青字部分がゼロになる）。よって、厳密解が存在するならば、の極限で、オイラー法はその解に収束する。