ニュートンのゆりかご（Newton's cradle）

 投稿者：山中和義  投稿日：2014年 7月13日(日)15時55分37秒
  おはじきのような球状の物体の運動と衝突を考えます。
シミュレーションには、ベクトル方程式を記述することにしました。
ベクトルの計算は、複素数平面すなわち複素数モードで計算ができますが、
数値がFPUによる浮動小数点のみになるので、このモードに限定しません。
十進モード、1000桁モード、2進モードに切替えて、実行してみます。

●その１　1次元（直線）の衝突
移動が、p=p+⊿pなので、誤差の累積が影響します。
十進モード、1000桁モードでは、うまく動作します。


!ニュートンのゆりかご（Newton's cradle）　○○○→　●●の場合

DIM WU(2),WD(2),WL(2),WR(2) !壁の表面の向き（法線ベクトル）
DATA 0,-1, 0,1, 1,0, -1,0
MAT READ WU !上
MAT READ WD !下
MAT READ WL !左
MAT READ WR !右

LET UP=10 !壁の位置（原点からの符号つき距離）　※
LET DW=-10
LET LF=-6
LET RT=14

SET WINDOW LF,RT,DW,UP !表示領域

PUBLIC NUMERIC E !跳ね返り係数
LET E=1 !弾性衝突のとき

!---------------------------------------

LET N=5 !ボールの個数　※
LET R=0.5 !ボールの半径　※

DIM X(N),Y(N) !位置(x,y)
DATA -4,0, -3,0, -2,0, 3,0, 4,0 !　※
FOR i=1 TO N
   READ X(i),Y(i)
NEXT i

DIM VX(N),VY(N) !移動速度vのX,Y成分
DATA 2,0, 2,0, 2,0, 0,0, 0,0 !　※
FOR i=1 TO N
   READ VX(i),VY(i)
NEXT i

!---------------------------------------

LET dt=0.1 !⊿t
DO
   SET DRAW mode hidden !ちらつき防止の開始
   CLEAR
   DRAW grid !座標

   FOR i=1 TO N !ボールを描く
      DRAW disk WITH SCALE(R)*SHIFT(X(i),Y(i))
   NEXT i

   SET DRAW mode explicit !ちらつき防止の終了

   !ボールの移動
   DO
      LET HIT=0

      FOR i=1 TO N !ボールiとボールjとの衝突を確認する
         DIM OA(2),AV(2), OB(2),BV(2)
         LET OA(1)=X(i) !copy it
         LET OA(2)=Y(i)
         LET AV(1)=VX(i)
         LET AV(2)=VY(i)

         IF DISTL(OA,WU,UP)<=R AND DOT(AV,WU)<0 THEN !表面から上壁に衝突するとき
            CALL HIT_WALL(AV,WU)
            LET HIT=-1
         END IF
         IF DISTL(OA,WD,-DW)<=R AND DOT(AV,WD)<0 THEN !下壁
            CALL HIT_WALL(AV,WD)
            LET HIT=-1
         END IF
         IF DISTL(OA,WL,-LF)<=R AND DOT(AV,WL)<0 THEN !左壁
            CALL HIT_WALL(AV,WL)
            LET HIT=-1
         END IF
         IF DISTL(OA,WR,RT)<=R AND DOT(AV,WR)<0 THEN !右壁
            CALL HIT_WALL(AV,WR)
            LET HIT=-1
         END IF

         FOR J=i+1 TO N !組み合わせによる
            DIM AB(2),Vd(2)
            LET OB(1)=X(J) !copy it
            LET OB(2)=Y(J)
            MAT AB=OB-OA !相対位置
            IF DOT(AB,AB)<=(2*R)^2 THEN !衝突するとき
               LET BV(1)=VX(J)
               LET BV(2)=VY(J)

               MAT Vd=AV-BV !ボールjを停止させた相対速度
               IF DOT(AB,Vd)>0 THEN !近づくとき（離れていくときは除く）

                  CALL HIT_BALL(OA,AV, OB,BV, R)

                  LET VX(J)=BV(1) !衝突後
                  LET VY(J)=BV(2)

                  LET HIT=-2
               END IF
            END IF
         NEXT J

         LET VX(i)=AV(1) !衝突後
         LET VY(i)=AV(2)
      NEXT i

      !!!PRINT HIT !debug
   LOOP UNTIL HIT=0 !衝突がなくなるまで

   FOR i=1 TO N
      LET X(i)=X(i)+dt*VX(i) !移動させる
      LET Y(i)=Y(i)+dt*VY(i)
   NEXT i

   WAIT DELAY 0.05 !アニメーションの速度  ※調整の必要性あり
LOOP
END

EXTERNAL SUB HIT_WALL(AV(), NW()) !ボールと壁との非弾性衝突による結果
DIM dV(2)
LET w=(1+E)*DOT(AV,NW) !壁の法線単位ベクトルnw
MAT dV=(w)*NW
MAT AV=AV-dV
END SUB

EXTERNAL SUB HIT_BALL(OA(),AV(), OB(),BV(), R) !同じ質量の2つのボールの非弾性衝突による結果
DIM OP(2),OC(2),V(2),dV(2)
MAT OP=OB-OA !衝突箇所の位置 ↑c
MAT OC=(1/2)*OP
MAT V=BV-AV !相対速度 v2-v1
LET w=(1+E)*DOT(OC,V)/(2*R*R) !w=(1+E)(↑c・(↑v2-↑v1))/(2r^2)
MAT dV=(w)*OC !ボールAの速度変化量
MAT AV=AV+dV !速度を変化させる
MAT BV=BV-dV
END SUB

EXTERNAL SUB VecNormalize(V()) !正規化する（長さが1の単位ベクトル）
LET L=SQR(DOT(V,V))
IF L<>0 THEN  MAT V=(1/L)*V
END SUB

EXTERNAL FUNCTION DISTL(P(),N(),D) !点pと直線↑n・↑P+d=0との距離
LET DISTL=ABS(DOT(N,P)+D)/SQR(DOT(N,N))
END FUNCTION



●その２　2次元（平面）の衝突

1000桁モード、2進モードでは、うまく動作します。

変更箇所

!---------------------------------------

LET N=2 !ボールの個数　※
LET R=1 !ボールの半径　※

DIM X(N),Y(N) !位置(x,y)
LET X(1)=-4
LET Y(1)=0
LET X(2)=SQR(3)
LET Y(2)=-1

DIM VX(N),VY(N) !移動速度vのX,Y成分
DATA 2,0, 0,0 !　※
FOR i=1 TO N
   READ VX(i),VY(i)
NEXT i

DEF F(X)=-1/SQR(3)*X !軌跡
DEF G(X)=SQR(3)*X

!---------------------------------------

LET dt=0.1 !⊿t
DO
   SET DRAW mode hidden !ちらつき防止の開始
   CLEAR
   DRAW grid !座標
   PLOT LINES: LF,F(LF); RT,F(RT)
   PLOT LINES: LF,G(LF); RT,G(RT)

   FOR i=1 TO N !ボールを描く
      DRAW disk WITH SCALE(R)*SHIFT(X(i),Y(i))
   NEXT i

   SET DRAW mode explicit !ちらつき防止の終了

   !ボールの移動

　　　以下省略
 

戻る