惑星の軌道

user-pic
0

もう30年以上前、20台後半でまだ独身の頃、シャープのMZという8ビット・パソコンを買い、BASICでプログラムを書いて、夜な夜な楽しく遊んでいたのだが、そのとき書いていたプログラムの中に、惑星の軌道計算があった。

惑星の軌道計算といっても、なんのことはない、ニュートンの法則の逆二乗則で力学の運動方程式を記述し、それを例によってルンゲ・クッタで数値計算するという素朴なものである。

それでも、惑星2つが引き合うという連星の二体問題を数値計算してプロットした結果を知り合いに見せたら、辛口で滅多に他人を褒めないその知り合いも、あまりこういうのは見たことがないので面白いと褒めてくれたので、うれしくなったのを覚えている。

それで気を良くして、ポアンカレによれば積分が存在しない、すなわち解析的に解けないといわれている三体問題をプロットしようとしたら、プログラムにバグがあるようでどうもうまくいかない。それでその後諦めていたのだが、最近、常微分方程式の数値解法に再び馴染んできたのをきっかけに、また再チャレンジしてみたくなった。

といっても、いきなり三体問題は荷が重いので、一体問題⇒二体問題⇒三体問題と進めていくことにした。一体問題というのはあまり聞きなれないが、私が意図したのは、中心に太陽があって、その周りを単一の惑星が周回する、というものである。しかしウェブで調べてみると、一体問題とは、宇宙にたった一つの星がある場合のことらしくて、するとニュートン力学のパラダイムでは、静止しているか等速直線運動するだけなので、数値計算もなにもあったものではない。

で、当初予定していたこのブログのタイトルを、一体問題から惑星の軌道という曖昧な言い方に変えて、中心に太陽があって、その周りを単一の惑星が周回する、という場合の数値計算をしてみることにした。

さて、それは要するに、太陽の周りに惑星が楕円軌道をとるだけだから楽勝と思っていたら、なかなか綾があった。すなわち、惑星を太陽の引力を受けて運動するようにするには、惑星の太陽からの距離と、太陽に対して相対的な速度を与えるのだが、ヨハン・ベルヌーイは、その場合惑星がとりえる軌道は、楕円、放物線、双曲線のどれかであることを見出した。

要するに、注意深く初期速度と位置を与えないと、軌道は放物線または双曲線となって、無限遠に「いっちゃう」のである。それで最初、適当に初期速度と位置与えて、軌道が発散して何度も困ってしまった。

さらに、自明なのでヨハン・ベルヌーイが指摘しなかったことがある。それはあまりに太陽の近くに初期位置を配置するか、初期速度が小さすぎると、惑星が太陽に引っ張られすぎて激突してしまうのである。これはこれでゲームセットである。

そうやってなんとか楕円軌道を得たのが、下記のプログラムで、(x,y,u,v) = (3.0, 2.0, 0.43. 0.01)の場合である。ここで、太陽を原点としたとき、xは惑星の初期x座標、yは惑星の初期y座標、uは初期速度のx方向成分、uは初期速度のy方向成分である。

その軌道は、以下のとおりだった。なお、下手なマンガみたいだが、太陽を赤い小円で示した。

これだと近日点が近すぎて、もし惑星が地球だと人類が焼け死んで滅びるかもしれない。そこで、妥当な距離の近日点になるように、初期位置と初期速度を、(x,y,u,v) = (3.0, 0.0, -0.1. 0.5)に変更してみた。すると軌道は、以下のとおりとなった。

このプログラムのソースコードを例によって以下に示す。それにしても、惑星や太陽の質量や重力定数やその他物理単位の次元を全く無視して適当に微分方程式を解くってどうよ、とつい自己突っ込みするのだった。

//---------------------------------------------------------------------------

#include <vcl.h>
#include <math.h>
#pragma hdrstop

#include "Main.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
double fx(double x, double y);
double fy(double x, double y);
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
	: TForm(Owner)
{
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Close1Click(TObject *Sender)
{

 Close();
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Run1Click(TObject *Sender)
{
const double H = 0.00001;  // 刻み
double x,y,u,v;     // x,yは座標値で、u,vはその速度のxとy成分
double a1,a2,a3,a4;  // dx
double b1,b2,b3,b4;  // dy
double c1,c2,c3,c4;  // du
double d1,d2,d3,d4;  // dv
double px, py;
int N = 0;           // 繰り返し回数

// 位置と速度の初期値
x = 3.0;  y = 0.0; u = -0.1; v = 0.5;

// 原点を描く
Canvas->Brush->Color = clRed;
Canvas->Ellipse(280,280,300,300);

while ( N++ < 10000000 )        {
 a1 = u;
 b1 = v;
 c1 = -fx(x,y);
 d1 = -fy(x,y);
 a2 = u + H*c1/2.0;
 b2 = v + H*d1/2.0;
 c2 = -fx(x+H*a1/2.0,y+H*b1/2.0);
 d2 = -fy(x+H*a1/2.0,y+H*b1/2.0);
 a3 = u + H*c2/2.0;
 b3 = v + H*d2/2.0;
 c3 = -fx(x+H*a2/2.0,y+H*b2/2.0);
 d3 = -fy(x+H*a2/2.0,y+H*b2/2.0);
 a4 = u + H*c3;
 b4 = v + H*d3;
 c4 = -fx(x+H*a3,y+H*b3);
 d4 = -fy(x+H*a3,y+H*b3);
 x += H*(a1 + 2.0*a2 + 2.0*a3 + a4)/6.0;
 y += H*(b1 + 2.0*b2 + 2.0*b3 + b4)/6.0;
 u += H*(c1 + 2.0*c2 + 2.0*c3 + c4)/6.0;
 v += H*(d1 + 2.0*d2 + 2.0*d3 + d4)/6.0;
 px = 290.0 + 70.0*x;
 py = 290.0 - 70.0*y;
 Canvas->Pixels[px][py] = clBlack;
}
}
//---------------------------------------------------------------------------
double fx(double x, double y) {
 double r,R;
 r = x*x + y*y;
 R = r * sqrt(r);
 return x/R;
}

double fy(double x, double y) {
 double r,R;
 r = x*x + y*y;
 R = r * sqrt(r);
 return y/R;
}

トラックバック(0)

トラックバックURL: http://commutative.world.coocan.jp/mt/mt-tb.cgi/4623

コメントする

このブログ記事について

このページは、あやたろうが2017年8月 8日 05:21に書いたブログ記事です。

ひとつ前のブログ記事は「ボーダーライン」です。

次のブログ記事は「きゅうりの船」です。

最近のコンテンツはインデックスページで見られます。過去に書かれたものはアーカイブのページで見られます。