微分方程式（数値解）
考えている物理的または化学的なシステムが複雑になると、 運動方程式などの微
分方程式で現象を記述できても、 その式を解析的に解くことはほとんど不可能に
なります。 そこで重要になってくるのが、微分方程式を数値的に解くテクニックです。
ここでは、 ２Dグラフィックスを使って式を解いた結果をビジュアルに示したいと思い
ます。勿論、三次元的な問題には３Dグラフィックスが必要になります。
まずは、仮想的に考えた太陽系を考えましょう。この太陽系には、恒星とその周りを
周回するいくつかの惑星が存在します。問題を簡単にするために、恒星の質量は惑
星の質量に比べて格段に大きいものとします（ 恒星の質量がそれほど大きくない場
合も後で検討します）。 この仮定によって、 惑星の軌道は惑星に関する運動方程
式のみを数値的に解けば、得られることになります。
例えば、 ある太陽系が一つの恒星と二つの惑星から構成されているとしましょう （
図１参照）。 このシステムでは、三つの運動方程式を同時に解かなければなりませ
ん。
恒星の運動方程式：
MS ｄ２RS /ｄｔ２＝−MS MP1 G（RS−RP１ ）/ |RS−RP１ |３
−MS MP２ G（RS−RP２ ）/ |RS−RP２ |３
惑星１の運動方程式：
MP１ ｄ２RP１ /ｄｔ２＝−MP１ MS G（RP１−RS ）/ |RP１−RS |３
−MP１ MP２ G（RP１−RP２ ）/ |RP１−RP２ |３
惑星２の運動方程式：
MP２ ｄ２RP２ /ｄｔ２＝−MP２ MS G（RP２−RS ）/ |RP２−RS |３
−MP２ MP１ G（RP２−RP１ ）/ |RP２−RP１ |３
ここで、MS、MP１、MP２、Gはそれぞれ、恒星の質量、惑星１の質量、惑星２の質量そ
して万有引力の定数です。
しかし、既に中心力の所で議論したように、MP１＜＜MSの場合は恒星の軌道の平
均半径が恒星自身の半径よりも格段に小さく、 惑星の円軌道（ または楕円軌道）の
中心（または焦点 ）は恒星の中心と一致していると考えて問題はありません。加えて
、もう一つの惑星による万有引力の影響も無視できるので、 解くべき微分方程式
は次の式だけになります。
惑星１の運動方程式：
MP１ ｄ２RP１ /ｄｔ２＝−MP１ MS G RP１ / |RP１ |３
上式におけるRP１ は恒星の中心から惑星１の中心に向かう位置ベクトルです。
では、上記の微分方程式を数値的に解くことを考えましょう。 まず、 その式を見れば
判るように、惑星１の加速度ベクトルは位置ベクトルに比例する形になっています。
したがって、 惑星１の速度ベクトルと位置ベクトルで決まる平面上に加速度ベクト
ルもあることが判ります。これは、惑星１が取る軌道が一つの平面上にあることを示
しています。 この平面上で微分方程式を考えれば、 問題自体は二次元の問題にな
ることに注意してください（図２参照）。
上で述べた平面上にXY軸をとり、運動方程式を成分表示してみます。
X成分： MP１ ｄ２XP１/ｄｔ２＝−MP１ MS G XP１/ （ X２P１＋Y２P１）３/２
Y成分： MP１ ｄ２YP１/ｄｔ２＝−MP１ MS G YP１/ （ X２P１＋Y２P１）３/２
この式において、座標の原点は恒星の中心と一致させてあります。円軌道の場合は
ともかく、XP１とYP1が二式の中でミックスしているため、簡単には解析的に解けませ
ん。以下に、極座標で表現した運動方程式も示しておきます。チャレンジ精神がある
方は自分で解いてみてください。これ以降は、惑星１の座標を表す添え字を省略しま
す。
微分方程式の挙動を調べる上で、係数などはあまり本質的ではないので質量や物
理定数は式から省くことにします。ゆえに、惑星１に関する運動方程式は、
X成分： ｄ２X/ｄｔ２＝−X/ （ X２＋Y２）３/２
Y成分： ｄ２Y/ｄｔ２＝−Y/ （ X２＋Y２）３/２
となります。以上から、解くべき微分方程式は上の二式になります。ただし、物理法則
を表す式では右辺と左辺の単位が揃わないといけないので、 気になる方は左辺に
単位を調節するための係数を導入しても構いません。
第一番目の式にYを掛けたものと第二番目の式にXを掛けたものは等しくなっている
ので、統合した形で微分方程式を表すと、
Y ｄ２X/ｄｔ２＝X ｄ２Y/ｄｔ２
となります。この式は、次のようにも書き換えられます。
（ｄ２X/ｄｔ２）/（ｄ２Y/ｄｔ２）＝X/Y
加速度ベクトルは位置ベクトルに比例していると指摘しましたが、上式はそのことをX
成分とY成分の比という形で示しています。
次に、数値計算するためのステップを考えてみます。初期位置（X０、Y０）と初期速度
（VX０、VY０）から、計算はスタートすることになります。時間的な刻み幅はΔｔとしまし
た。
ステップ１
（X０、Y０）から、運動方程式を使って初期加速度（AX０、AY０）を計算する。
AX０＝−X０/ （ X０２＋Y０２）３/２、 AY０＝−Y０/ （ X０２＋Y０２）３/２
ステップ２
（X０、Y０）と（VX０、VY０）から、Δｔ後の（X１、Y１）を計算する。
X１＝X０＋VX０Δｔ、 Y１＝Y０＋VY０Δｔ
ステップ３
（VX０、VY０）と（AX０、AY０）から、Δｔ後の（VX1、VY１）を計算する。
VX１＝VX０＋AX０Δｔ、 VY１＝VY０＋AY０Δｔ
後は、このプロセスを繰り返すことによって、 （X２、Y２）、（X３、Y３）、（X４、Y４）・・・・・・
という具合に順番に座標を計算して行けば良い訳です。 これらの点を２Dグラフィック
スを使って全て繋いでやれば、最終的に惑星１の軌道になります。
それでは、初期値に数値を与えて実際に軌道を計算してみます。 使うソフトは十進B
ASICです。以下にそのプログラムを示します。
上のプログラムにおいて、 初期値を（X０、Y０）＝（２０、０）、（VX０、VY０）＝（０、０．２
８）とし、時間ステップを３０００まで取った結果を以下に示します。
また、上の数値を２Dグラフックスを使って表示したものは次のようになります。
計算の精度に問題があるため軌道は閉じていませんが（ 図３参照、変更後 ）、プ
ログラムを工夫すればより正確な軌道を計算して描画させることは可能です。
私たちが住む太陽系に関しては、太陽や地球の質量は判っています。それらを基に
、初期位置や初期速度をいろいろ仮定して地球の軌道を数値計算してみてください
。