この実習では Processing を使用します。 Processing は、高度なプログラム言語であるＣやＣ＋＋と同様の言語体系を持ち、 同時に容易にグラフィックを描ける利便性を兼ね備えており、 プログラムの初歩を学ぶにはとても適切なアプリケーションです。 簡単な目的の計算なら十分にこなす事ができ、 また高度なプログラム環境に移行する必要が生じても、 Processing で勉強した事が無駄になりません。 皆さんの今後の学習にきっと役立つと思います。
この実習では、Processing の使い方の初歩から始めて、 物理学の基本である、ニュートン方程式を数値的に解き、 色々な物理現象が１つの基本原理と１つの解法で説明できる事を体験して頂きます。
情報教育システムにはV3.3.5のProcessingが既にインストールされています。 この実習ではこれを利用して頂きます。
もし自宅のＰＣ等で使用したい場合は、オリジナルサイトから ダウンロードして使う事も可能です。 但し、オリジナルサイトでは最新のV3.3.7があります。バージョンが異なると動作や使い方が変わることがあるので、注意が必要です。
また、レポートで使用するバージョンは情報教育システムと同じバージョンとします。 これは公平性を保つためです。 皆さんが各自でより進んだ内容を勉強する事は大変よいことではありますが、 レポートで使用できるのはV3.3.5である事を忘れないで下さい。
クイック起動バーから「Devel」→「processing」を選択して起動します。
中ほどの白い部分にプログラムを書く。 まずは以下のプログラムを打ち込んでみよ。構文に慣れるために、コピー＆ペーストは使用しないこと。
int i;
for(i=0; i<=10; i++){
println(i);
}
左上に、丸の中に三角があるボタンが、プログラムを走らせるボタンである。押して見よ。
下方の黒い部分に文字の出力が行われる。ここに、０から１０の整数が表示された事を確認せよ。
先程のボタンの１つ右にある、丸の中に四角があるボタンがプログラムを止めるボタンである。 押して、プログラムを止めよ。
プログラムの中にある for 文は、繰り返しを行うための構文である。 括弧の中は、左から、
もし終了条件を絶対に満たさないようなプログラムにして実行してしまうと、 当然であるが、プログラムは止まらない。この状況になっても「停止ボタン」を 押して、止める事ができる。
入力している途中で、灰色の部分が赤くなり、そこにエラーメッセージが表示される場合がある。 人間が入力している途中では、一時的に間違った文になるのは当たり前なので、気にしなくてよい。 但し、入力が完了した場合でもエラーメッセージが出ている場合は、それをよく読んで対応すること。 しかし、エラーメッセージが常に適切な指摘をしてくれる訳ではない事にも留意せよ。
尚、初心者はエラーメッセージの内容が理解できない事が多い。 恐らくこういう意味だろうと思って試してみるのは良い事であるが、 一方で、やみくもに多くを試してみるのは効率的ではない。 意味をよく推察してみて、それでも分からない場合は、ＴＡに聞いて下さい。
以下のプログラムを入力して実行してみよ。 コピー＆ペーストは使用しないこと。
int i;
size(800, 600);
for(i=0; i<350; i++){
println("X=", i*2, "Y=", i);
point( i*2, i);
}
size 文で、描画する範囲を横８００ピクセル、縦６００ピクセルに指定している。
i を０から３４９まで、１つづ増やしながら、 Ｘ座標はｉの２倍、Ｙ座標はｉの位置に、点を描画している。 Ｘ，Ｙの値は数値でも出力している。
点が小さくて見にくい場合は、point の代わりに、ellipse を使うと見やすくなる。
ellipse( i*2, i, 1, 1);
ellipse は楕円を描画する関数である。半径の小さな楕円を描く事で、大きな点を描画している。
以下の事をやってみよ。
尚、Processingで使用できる関数等の説明が
以下のプログラムを入力して実行してみよ。 コピー＆ペーストは使用しないこと。
void setup() {
size(600, 600);
println("width=", width, "height=", height);
frameRate(30);
line(sx(-100), sy(0), sx(+100), sy(0));
line(sx(0), sy(-100), sx(0), sy(+100));
}
float x1=-100.0, y1=-100.0;
float x2=+100.0, y2=+100.0;
float sx(float x){ return width *(x-x1)/(x2-x1); }
float sy(float y){ return height*(y2-y)/(y2-y1); }
float x=0.0, y;
void draw(){
y=0.01*x*x;
println("X=", x, "Y=", y);
point( sx(x), sy(y)); // ellipse(sx(x),sy(y),1,1);
x=x+1.0;
if(x>100.0) noLoop();
}
９～１３行は、座標変換のお膳立てをしている。 size 文で宣言された６００×６００ピクセルのスクリーンは、いわばＸが１～６００、Ｙが１～６００の 座標を持っていると言える。但し、Ｙは下が大きい事に留意せよ。 これだと我々が自由に座標系を使えない。そこで、計算で使用する座標系からスクリーンに表示する座標系への 変換を考える。このプログラムでは計算で使用するＸは、－１００．０から＋１００．０まで、 Ｙも同じ範囲とし、Ｙ軸は上が大きくなるように方向を逆転させる、一次変換を行う関数を定義している。
この関数はＸ，Ｙ毎に１つづつ定義されている。 計算に使用する座標系は９～１０行目で定義されている。 変換に必要なスクリーンのピクセル数は、size 文の実行後に、自動的に、width, height という変数に代入されている。
この関数は、下から４行目の point 関数で使用されている。 この様にすると、point関数が計算座標系のＸ、Ｙを引数に持っているかのごとく使える。
プログラムでは、２次関数を計算、描画している。
Ｘを０から１づつ増やし、１００を超えたらループを終了している。 Ｘの初期値、終了条件、増分、関数形等を変更して、色々試してみよ。
setup関数はプログラムの最初に１度だけ実行される。
draw関数は、何もしなければ、自動的に無限回、実行される。 繰り返し計算し、その結果を描画する場合には、draw関数の中に、１回分の描画するプログラムを書けば、 それを繰り返してくれる。但し、パラメータの更新や、終了条件の判断は自分で書く必要がある。
関数形を変えて、異なるグラフを描け。
初期値、終わり値を変えてみよ。
スクリーンの大きさを変えてみよ。また、計算に使う座標系も変えてみよ。
draw関数はsetup関数の中でframeRate関数で設定された頻度で繰り返される。 このプログラムでは３回／秒としてある。この数値を変えてみよ。
以下のプログラムを入力して実行してみよ。 これ以降、コピー＆ペーストを使用してもよいが、 プログラムの中身をしっかりと 読んで、理解する事。
void setup() {
size(600, 600);
println("width=", width, "height=", height);
//frameRate(3);
line(sx(-1000), sy(0), sx(+1000), sy(0));
line(sx(0), sy(-1000), sx(0), sy(+1000));
}
float x1= -1.0, y1=-100.0;
float x2=+10.0, y2= +10.0;
float sx(float x){ return width *(x-x1)/(x2-x1); }
float sy(float y){ return height*(y2-y)/(y2-y1); }
int i=0;
float t=0.0, dt=0.01;
float x=0.0, y=0.0;
float vx=10.0, vy=10.0;
float g=-9.8;
void draw(){
t=t+dt;
vx=vx;
vy=vy+g*dt;
x=x+vx*dt;
y=y+vy*dt;
println("i=", i, "t=", t, "Vx=", vx, "Vy=", vy, "X=", x, "Y=", y);
//point( sx(t), sy(vx));
//point( sx(t), sy(vy));
//point( sx(t), sy(x));
point( sx(t), sy(y));
//point( sx(x), sy(y));
if(i++>500) noLoop();
}
このプログラムは、質点の放物運動を数値計算で解いている。
物理学で言う、数値計算で問題が「解けた」とは、
重力加速度は、－９．８ｍ／ｓ／ｓ（下方へなので符号は負）とし、 質点を座標（０ｍ，０ｍ）から初速度（１０ｍ／ｓ、１０ｍ／ｓ）で飛ばした時の運動を追跡している。
時間の刻みは０．０１秒とし、ある時刻から次の時刻の速度と座標は以下のように求めている。
Vx ← Vx
Vy ← Vy + g・Δt
X ← X + Vx・Δt
Y ← Y + Vy・Δt
これらの式は、速度と加速度の定義から分かると思う。
Vx(t) = ( X(t+Δt) - X(t) ) / Δt （Δt→０）
Ax(t) = ( Vx(t+Δt) - Vx(t) ) / Δt （Δt→０）
この様に１次の極限での近似を用いて数値積分する方法を、オイラー法という。
プログラムの最後の方に、Ｘ軸とＹ軸に何を描画するかが、何行かに渡り書いてあり、 １つを除きコメントアウトされている。これらをそれぞれ描画し、計算の結果を確認せよ。 初速度を変えて試してみよ。
このように、物理現象をシミュレーションする場合には、物理法則に関わる部分には一切修正を加えず、 初期値等の条件を変更する事で、様々な物理現象を再現できる。
Δtは０．０１秒で計算したが、これを０．０２秒に粗くしてみよ。 グラフでは分からない程度の差しかでないので、出力された数値を比較して、計算誤差が変化する事を確認せよ。
// （スラッシュ２個）は、コメント機能である。この記号から行末までは、コメントとなり、プログラムの実行に何の影響も与えない。 また、/* コメント */ という構造を使用すれば、複数行にまたがったコメントも作成できる。
以下のプログラムを入力して実行してみよ。
void setup() {
size(600, 600);
println("width=", width, "height=", height);
//frameRate(1000);
colorMode(RGB);
background(0,0,0);
stroke(255,0,0);
fill(255,0,0);
ellipse( sx(0), sy(0), 5, 5);
stroke(0,255,0);
fill(0,255,0);
}
float x1=-1.6e11, y1=-1.6e11;
float x2=+1.6e11, y2=+1.6e11;
float sx(float x){ return width *(x-x1)/(x2-x1); }
float sy(float y){ return height*(y2-y)/(y2-y1); }
int i=0;
float t=0.0, dt=80000.0;
float x=1.496e11, y=0.0;
float vx=0.0, vy=2.979e4;
float G=6.67384e-11;
float Ms=1.989e30;
void draw(){
float r3;
t=t+dt;
r3=pow( sq(x)+sq(y), 1.5);
vx=vx-dt*G*Ms*x/r3;
vy=vy-dt*G*Ms*y/r3;
x=x+vx*dt;
y=y+vy*dt;
println("i=", i, "t=", t, "Vx=", vx, "Vy=", vy, "X=", x, "Y=", y);
point( sx(x), sy(y));
if(i++>10000) noLoop();
}
このプログラムは、太陽の周りを回る地球の動きをシミュレーションしている。
最初の方で、色を指定している。
初期条件として、太陽と地球の距離、地球の速度は、現実的な値を入れてあるが、 この場合、綺麗に円軌道が再現できる。
初速度を変えると楕円軌道になる筈なので、やってみよ。
初速を半分（1.5e4）にすると楕円軌道がだんだんずれていってしまう。 これは計算誤差が蓄積するからである。ｄｔの値を小さくすると計算誤差が減り、 より正しい結果を得られる事を確かめよ。。
ｄｔの値を小さくした場合は、frameRate で大きな値を指定すると速く計算する。 指定しなかった時の値は６０である。
以下のプログラムを入力して実行してみよ。 コピー＆ペーストを使用してもよい。
void setup() {
size(600, 600);
println("width=", width, "height=", height);
//frameRate(1000);
background(0);
stroke(255);
colorMode(HSB,360,100,100);
}
final float G=6.674e-11;
/////////// Parameters to change an initial condition of the probe ///////////
float theta=PI/3; // Direction of momentum.
float disp=1.8e10; // Radial displacement between the Jupiter and the probe.
//////////////////////////////////////////////////////////////////////////////
float x1=+6.784e11, y1=-1.e11;
float x2=+8.784e11, y2=+1.e11;
/*
float x1=-2.e12, y1=-2.e12;
float x2=+2.e12, y2=+2.e12;
*/
float sx(float x) {
return width *(x-x1)/(x2-x1);
}
float sy(float y) {
return height*(y2-y)/(y2-y1);
}
float accel(float m, float r) {
return G*m/sq(r);
}
int i=0;
float t=0.0, dt=10000.0;
float x_jup=7.784e11*cos(-PI/30), y_jup=7.784e11*sin(-PI/30);
float x_prb=7.784e11-disp, y_prb=-7.e10;
float vx_jup=1.306e4*(-sin(-PI/30)), vy_jup=1.306e4*cos(-PI/30);
float vx_prb=1.e4*cos(theta), vy_prb=1.e4*sin(theta);
float m_sol=1.989e30;
float m_jup=1.898e27;
float m_prb=721.9; //voyger 1 & 2
float ax_jup=0.0, ay_jup=0.0;
float ax_prb=0.0, ay_prb=0.0;
// the sun stands still //
void draw() {
float r_sol_jup, r_sol_prb, r_jup_prb;
t=t+dt;
r_sol_jup=sqrt(sq(x_jup)+sq(y_jup));
r_sol_prb=sqrt(sq(x_prb)+sq(y_prb));
r_jup_prb=sqrt(sq(x_jup-x_prb)+sq(y_jup-y_prb));
ax_jup=-x_jup/r_sol_jup*accel(m_sol, r_sol_jup);
ay_jup=-y_jup/r_sol_jup*accel(m_sol, r_sol_jup);
ax_prb=(x_jup-x_prb)/r_jup_prb*accel(m_jup, r_jup_prb)-x_prb/r_sol_prb*accel(m_sol, r_sol_prb);
ay_prb=(y_jup-y_prb)/r_jup_prb*accel(m_jup, r_jup_prb)-y_prb/r_sol_prb*accel(m_sol, r_sol_prb);
vx_jup=vx_jup+dt*ax_jup;
vy_jup=vy_jup+dt*ay_jup;
vx_prb=vx_prb+dt*ax_prb;
vy_prb=vy_prb+dt*ay_prb;
x_jup=x_jup+dt*vx_jup;
y_jup=y_jup+dt*vy_jup;
x_prb=x_prb+dt*vx_prb;
y_prb=y_prb+dt*vy_prb;
stroke(0,0,99);
point( sx(x_jup), sy(y_jup));
stroke(180,99,99);
point( sx(x_prb), sy(y_prb));
println("V=", sqrt(sq(vx_prb)+sq(vy_prb)), "ALT=", r_jup_prb);
if (i++>10000) noLoop();
}
このプログラムは、スイングバイをシミュレーションしている。 スイングバイについては、「講義資料（物理）追加」を参照せよ。
課題５は２体問題（２つの質点の運動を扱う）であったが、これは、太陽、木星、探査機の ３体問題である。尚、数学的には３体以上の厳密解はなく、高精度の結果を得るには数値的に 解くしか方法がない。また、３体ではあるが、太陽は座標原点に固定してある事に留意せよ。
「講義資料（物理）追加」を参照して、初期条件を変え、スイングバイが色々なパターンを 取る事を確認せよ。地球に帰還する軌道に載せたり、惑星軌道に載せたり、 次の惑星に向かったり、太陽系を脱出したりする事ができる。
以下のプログラムを入力して実行してみよ。 コピー＆ペーストを使用してもよい。
void setup() {
size(600, 600);
println("width=", width, "height=", height);
frameRate(1000);
colorMode(RGB);
background(0,0,0);
stroke(255,0,0);
fill(255,0,0);
ellipse( sx(0), sy(0), 5, 5);
stroke(0,255,0);
fill(0,255,0);
}
float x1=-1.6e11, y1=-1.6e11;
float x2=+1.6e11, y2=+1.6e11;
float sx(float x){ return width *(x-x1)/(x2-x1); }
float sy(float y){ return height*(y2-y)/(y2-y1); }
int i=0;
float t=0.0, dt=100.0; // dt=10.0 or 1.0
float x=1.496e11, y=0.0;
float vx=0.0, vy=3.0e4;
float G=6.67384e-11;
float Ms=1.989e30;
void draw(){
float r, r3;
for(int j=0; j<10000; j++){
t=t+dt;
r =sqrt(x*x+y*y);
r3=r*r*r;
vx=vx-dt*G*Ms*x/r3;
vy=vy-dt*G*Ms*y/r3;
x=(x+vx*dt);
y=(y+vy*dt);
}
println("i=", i, "t=", t, "Vx=", vx, "Vy=", vy, "X=", x, "Y=", y);
point( sx(x), sy(y));
if(i++>10000) noLoop();
}
このプログラムは課題５を元にして修正を加えてある。
より正確な計算をする目的で、ｄｔの値を１００に小さくしてある。 普通に計算すると８００倍の計算時間がかかってしまう。 しかし、計算時間の大部分は画面表示にかかっている。 １回の繰り返しで１回の画面表示を行うと時間がかかるので、１００００回に１回の画面表示を行うように 修正してある。積分には実際に８００倍の計算を行わせている事に留意せよ。
このまま実行させ、円軌道が再現されている事を確認せよ。 また実行時間を把握せよ。
次に、ｄｔを１０にして実行してみよ。 同じ軌道を再現するのに１０倍の時間がかかり、より高精度で計算ができた筈である。
次に、ｄｔを１にして実行してみよ。 同じ軌道を再現するのに更に１０倍の時間がかかり、更に高精度で計算ができた筈であるが、、、
軌道の最初の部分が直線になってしまっている。何が起きているのか、調べるためにプログラムを再度起動し、 すぐに（１秒くらい）止めよ。コンソールをスクロールアップし、最初の部分を見て、Ｘの値が変化していない 事に気が付け。それは何故か考えよ。
以下のプログラムを入力して実行してみよ。 コピー＆ペーストを使用してもよい。
void setup() {
size(600, 600);
println("width=", width, "height=", height);
frameRate(1000);
colorMode(RGB);
background(0,0,0);
stroke(255,0,0);
fill(255,0,0);
ellipse( sx(0), sy(0), 5, 5);
stroke(0,255,0);
fill(0,255,0);
}
float x1=-1.6e11, y1=-1.6e11;
float x2=+1.6e11, y2=+1.6e11;
float sx(float x){ return width *(x-x1)/(x2-x1); }
float sy(float y){ return height*(y2-y)/(y2-y1); }
double dsqrt(double x){
double x0, x1;
x0=x; x1=(x+1.0D)/2.0D;
while (x1!=x0) {
x0=x1;
x1=(x/x0+x0)/2.0D;
}
return x1;
}
int i=0;
double t=0.0, dt=1.0; // dt=1.0
double x=1.496e11, y=0.0;
double vx=0.0, vy=3.0e4;
double G=6.67384e-11;
double Ms=1.989e30;
void draw(){
double r, r3;
for(int j=0; j<10000; j++){
t=t+dt;
r =dsqrt(x*x+y*y);
r3=r*r*r;
vx=vx-dt*G*Ms*x/r3;
vy=vy-dt*G*Ms*y/r3;
x=(x+vx*dt);
y=(y+vy*dt);
}
println("i=", i, "t=", t, "Vx=", vx, "Vy=", vy, "X=", x, "Y=", y);
point( sx((float)x), sy((float)y));
if(i++>10000) noLoop();
}
このプログラムは、課題７のプログラムを元にして、倍精度で計算するように変更してある。 下記のようにプログラムの核心の部分を変更せず、単に計算精度を上げただけで、問題なく計算ができるようになった事を確認せよ。
３０～３４行目と３７行目では、計算に使用する変数を、単精度（float）から倍精度（double）に変更して宣言してある。 このように宣言を変更するだけで、実際に計算するところの変更は必要がない。
描画に使用する部分は単精度のままで変更していない（１４～１７行目）。 但し、５０行目では、単精度でしか受け付けない関数 sx, sy に引き渡す際に、倍精度から単精度に変換している。
平方根を求める関数 sqrt は、倍精度で計算する事ができないので、 １９～２７行目で、倍精度で平方根を計算する関数 dsqrt を定義している。 数字の後ろに付いている D は、それが倍精度の数値である事を示している。
この実習では、単に倍精度の平方根を求める関数が定義され使われている事を理解すればよいが、 余力のあるものは、このアルゴリズムで平方根が求められる事を証明してみるのもよいだろう。 更に余力のあるものは、このプログラムでは平方根が求められない場合がある事を示し、 その欠点を改良したプログラムを考えてみるのも良いだろう。
倍精度を使わずに、単精度のままであっても、アルゴリズムを改良する事で計算精度を上げる方法もある。
まず、課題５のプログラムを修正して、ｄｔを１００００、ｖｙの初期値を、１．０ｅ４にして、 実行してみよ。実行に時間がかかるようになるので、コメントアウトしてある frameRate(1000) を有効にして、 速く計算させよ。 楕円軌道がだんだんずれていく事が分かる。厳密な解ではこれは起きない。
次に、同じ初期条件と時間ステップで、４次のルンゲクッタ法を用いたプログラムを実行してみる。 以下のプログラムを入力して実行してみよ。 コピー＆ペーストを使用してもよい。
void setup() {
size(600, 600);
println("width=", width, "height=", height);
frameRate(1000);
colorMode(RGB);
background(0,0,0);
stroke(255,0,0);
fill(255,0,0);
ellipse( sx(0), sy(0), 5, 5);
stroke(0,255,0);
fill(0,255,0);
}
float x1=-1.6e11, y1=-1.6e11;
float x2=+1.6e11, y2=+1.6e11;
float sx(float x){ return width *(x-x1)/(x2-x1); }
float sy(float y){ return height*(y2-y)/(y2-y1); }
int i=0,j=0;
float[][] k = new float[4][4];
float t=0.0, dt=10000.0;
float x=1.496e11, y=0.0;
float vx=0.0, vy=1.e4;
float ax=0.0, ay=0.0;
float G=6.67384e-11;
float Ms=1.989e30;
void draw(){
float r3;
t=t+dt;
//step 1
r3=pow( sq(x)+sq(y), 1.5);
ax=-G*Ms*x/r3;
ay=-G*Ms*y/r3;
k[0][0]=vx*dt;
k[2][0]=ax*dt;
k[1][0]=vy*dt;
k[3][0]=ay*dt;
//step 2
r3=pow( sq(x+0.5*k[0][0])+sq(y+0.5*k[1][0]), 1.5);
ax=-G*Ms*(x+0.5*k[0][0])/r3;
ay=-G*Ms*(y+0.5*k[1][0])/r3;
k[0][1]=(vx+0.5*k[2][0])*dt;
k[2][1]=ax*dt;
k[1][1]=(vy+0.5*k[3][0])*dt;
k[3][1]=ay*dt;
//step 3
r3=pow( sq(x+0.5*k[0][1])+sq(y+0.5*k[1][1]), 1.5);
ax=-G*Ms*(x+0.5*k[0][1])/r3;
ay=-G*Ms*(y+0.5*k[1][1])/r3;
k[0][2]=(vx+0.5*k[2][1])*dt;
k[2][2]=ax*dt;
k[1][2]=(vy+0.5*k[3][1])*dt;
k[3][2]=ay*dt;
//step 4
r3=pow( sq(x+k[0][2])+sq(y+k[1][2]), 1.5);
ax=-G*Ms*(x+k[0][2])/r3;
ay=-G*Ms*(y+k[1][2])/r3;
k[0][3]=(vx+k[2][2])*dt;
k[2][3]=ax*dt;
k[1][3]=(vy+k[3][2])*dt;
k[3][3]=ay*dt;
//last
x=x+(k[0][0]+2*k[0][1]+2*k[0][2]+k[0][3])/6;
y=y+(k[1][0]+2*k[1][1]+2*k[1][2]+k[1][3])/6;
vx=vx+(k[2][0]+2*k[2][1]+2*k[2][2]+k[2][3])/6;
vy=vy+(k[3][0]+2*k[3][1]+2*k[3][2]+k[3][3])/6;
println("i=", i, "t=", t, "Vx=", vx, "Vy=", vy, "X=", x, "Y=", y);
point( sx(x), sy(y));
if(i++>10000) noLoop();
}
ルンゲクッタ法の場合は、見える範囲内では楕円はずれず、計算精度が上がっている事が分かる。
一方で計算量は大幅に増えているが、画面表示にかかる時間の方が大きいので、実際の待ち時間は 変わらない。
オイラー法でもｄｔを２０００とすれば楕円がずれない精度が出せる。 しかし同じ時間をシミュレートするための計算量は５倍となる。
オイラー法は、計算量が少なく、精度が低い。ルンゲクッタ法は、計算量が多く、精度が高い。 今回の例では非常に簡単な例であったため、両者の違いは大きくなかったが、 ほとんどの科学技術計算において、ルンゲクッタ法（更にはその改良版）が使われている。
以下のプログラムを入力して実行してみよ。 コピー＆ペーストを使用してもよい。
int i=0,j;
float t=0.0, dt=10;
float x[]={0,0,0}, y[]={0,0,0};
float vx[]={0,0,0}, vy[]={0,0,0};
float ax[]={0,0,0}, ay[]={0,0,0};
float r[][]={{0,0,0},{0,0,0},{0,0,0}};
float m[]={5.e30,5.e30,5.e30};
void setup() {
x[0]=0.;
y[0]=5.e10;
x[1]=x[0]*cos(2*PI/3)-y[0]*sin(2*PI/3);
y[1]=x[0]*sin(2*PI/3)+y[0]*cos(2*PI/3);
x[2]=x[0]*cos(-2*PI/3)-y[0]*sin(-2*PI/3);
y[2]=x[0]*sin(-2*PI/3)+y[0]*cos(-2*PI/3);
vx[0]=5.e4;
vy[0]=2.e4;
vx[1]=vx[0]*cos(2*PI/3)-vy[0]*sin(2*PI/3);
vy[1]=vx[0]*sin(2*PI/3)+vy[0]*cos(2*PI/3);
vx[2]=vx[0]*cos(-2*PI/3)-vy[0]*sin(-2*PI/3);
vy[2]=vx[0]*sin(-2*PI/3)+vy[0]*cos(-2*PI/3);
size(600, 600);
println("width=", width, "height=", height);
background(0);
stroke(255);
colorMode(HSB,360,100,100);
}
final float G=6.674e-11;
float x1=-1.e11, y1=-1.e11;
float x2=+1.e11, y2=+1.e11;
float sx(float x) {
return width *(x-x1)/(x2-x1);
}
float sy(float y) {
return height*(y2-y)/(y2-y1);
}
float accel(float m, float r) {
return G*m/sq(r);
}
void draw() {
for(j=0;j<500;j++){
t=t+dt;
r[0][1]=sqrt(sq(x[0]-x[1])+sq(y[0]-y[1]));
r[1][2]=sqrt(sq(x[1]-x[2])+sq(y[1]-y[2]));
r[2][0]=sqrt(sq(x[2]-x[0])+sq(y[2]-y[0]));
ax[0]=(x[1]-x[0])/r[0][1]*accel(m[1],r[0][1])+(x[2]-x[0])/r[2][0]*accel(m[2],r[2][0]);
ay[0]=(y[1]-y[0])/r[0][1]*accel(m[1],r[0][1])+(y[2]-y[0])/r[2][0]*accel(m[2],r[2][0]);
ax[1]=(x[2]-x[1])/r[1][2]*accel(m[2],r[1][2])+(x[0]-x[1])/r[0][1]*accel(m[0],r[0][1]);
ay[1]=(y[2]-y[1])/r[1][2]*accel(m[2],r[1][2])+(y[0]-y[1])/r[0][1]*accel(m[0],r[0][1]);
ax[2]=(x[0]-x[2])/r[2][0]*accel(m[0],r[2][0])+(x[1]-x[2])/r[1][2]*accel(m[1],r[1][2]);
ay[2]=(y[0]-y[2])/r[2][0]*accel(m[0],r[2][0])+(y[1]-y[2])/r[1][2]*accel(m[1],r[1][2]);
vx[0]=vx[0]+dt*ax[0];
vy[0]=vy[0]+dt*ay[0];
vx[1]=vx[1]+dt*ax[1];
vy[1]=vy[1]+dt*ay[1];
vx[2]=vx[2]+dt*ax[2];
vy[2]=vy[2]+dt*ay[2];
x[0]=x[0]+dt*vx[0];
y[0]=y[0]+dt*vy[0];
x[1]=x[1]+dt*vx[1];
y[1]=y[1]+dt*vy[1];
x[2]=x[2]+dt*vx[2];
y[2]=y[2]+dt*vy[2];
}
stroke(0,99,99);
point( sx(x[0]), sy(y[0]));
stroke(120,99,99);
point( sx(x[1]), sy(y[1]));
stroke(240,99,99);
point( sx(x[2]), sy(y[2]));
println("X0=", x[0], "Y0=", y[0], "X1=", x[1], "Y1=", y[1], "X2=", x[2], "Y2=", y[2]);
if (i++>10000) noLoop();
}
このプログラムは、連星系のシミュレーションをしている。
課題６では太陽の座標を原点に固定していたが、この計算では３体とも固定されておらず、 全て計算している。
配列が使用されている。
高精度の計算が必要になってくるので、ｄｔは非常に小さくなっている。 そのため、５００回の繰り返しを計算して、１回の描画を行うようになっている。
ニュートンの万有引力の法則下では、２体問題（解くべき物体が２個）の場合は、数学的な厳密解があり、 それが楕円軌道（と彗星のような双曲線軌道）になる事が分かっている。しかし、３体以上になると 数学的な厳密解がない事も分かっている。（解がないという意味ではない。式に書けないという意味である。） このような場合、数学的に近似する（つまりモデルを変更する）方法もあるが、 それを行ってしまうと得られた数値の精度がモデルの精度以上にはできないという限界がある。 より高精度な結果を得るには、数学的な近似をせず、計算上の近似を高めていく手法しかない。 それがシミュレーションが必要な理由の１つである。
ＣＬＥにレポート問題を出してある。 問題とその注意をよく読んで、提出せよ。 レポートは単位の認定に必要なものである。