gravity simulation #1

シミュレーション系の手始めに二次元で．
質点を一個原点に固定して，もう一個の質点を初速度与えてその後の運動をシミュレート．

function grav2d;
nstep=10000;  //step 数
dt=0.01;//時間刻み
r=zeros(nstep+1,2);  //位置
v=zeros(nstep+1,2);  //速度
a=zeros(nstep+1,2);  //加速度
GM=input('gm?');  //固定した質点のM と重力定数の積
r(1,:)=input('r1?');  //初期値設定
v(1,:)=input('v1?');
a(1,:)=r(1,:)*GM/(norm(r(1,:))^3);
for n=2:nstep+1,
    r(n,:)=r(n-1,:)+dt*v(n-1,:);
    v(n,:)=v(n-1,:)+dt*a(n-1,:);
    a(n,:)=-r(n,:)*GM/(norm(r(n,:))^3);
    if norm(r(n,:))<0.1 then break;
    end
end
x=r(:,1); y=r(:,2);   //こことその下は一行で書けるのかも
plot2d(x,y);

norm の三乗とかでるからどうしても ちょっと計算量は多くなる．
なおまだちゃんと数値を考えた上で動かして見てないから正しく動いてるかちょっと不明．

<追記 @15:21 01/12/2011>
このプログラムのままでは，例えば円軌道を描くべき初期設定でも楕円を描いてしまう．
dt が大きすぎるため，と思われ，10^{-4} くらいでちゃんと動くようだ(thanks; @sericum )．
ついでに，scilab の Floating-point relative accuracy は %eps .