2重振子メーター付(解析したい方へ)再投稿

 投稿者:SECOND  投稿日:2009年 1月10日(土)02時54分20秒
  !2重振子メーター付(解析したい方へ)再投稿

!2重振り子特有の下降衝撃時に、その急変する過程が、
!演算ピッチの間に入って脱落し、全エネルギーが変動する計算エラーが見られる。
!演算ピッチが小さいと、緩和するが、低速パソコンでは、描画ピッチが伴わない。

!-----
LET g= 9.8 ! m/s^2 重力加速度
LET m1=.188 ! kg おもり
LET m2=.188 ! kg
LET L1= 5 ! m 吊り棒
LET L2= 5 ! m
LET r1=.75*SQR(m1) ! おもりの描画径
LET r2=.75*SQR(m2)
!
LET dt=0.05 !sec. 演算ピッチ。高速機 ほど、小さく。(0.05は、Pentium3 500MHz)
!
!※0.01くらいが望ましいが、描画ピッチ(画面に表示)が、ついて来れなくなったら戻す。
! 2つのピッチが、ズレていると物理的な速度ではなくなる。
!
!------------ 2重振り子の方程式
! d(dθ1)/dt^2=
! [ g*{sinθ2*cosΔ-μ*sinθ1}-{L2*(θ2/dt)^2+L1*(θ1/dt)^2*cosΔ}*sinΔ]
!          /[ L1*{μ-cosΔ^2}]
! d(dθ2)/dt^2=
! [ g*μ*{sinθ1*cosΔ-sinθ2}+{μ*L1*(θ1/dt)^2+L2*(θ2/dt)^2*cosΔ}*sinΔ]
!          /[ L2*{μ-cosΔ^2}]
!
! g= 重力加速度 μ=(m1+m2)/m2 Δ=θ1-θ2
!
!---------- 式の整理(θ1θ2 共に、重力方向0からの左回り角)
LET μ2=m2/(m1+m2)
LET L21=L2/L1
!
!ss1=-g/L1*sinθ1 -μ2*L21*ω2^2*sinΔ
!ss2=-g/L2*sinθ2 +    ω1^2*sinΔ/L21
!D=1-μ2*COSΔ^2
!d(θ1)/dt=ω1
!d(θ2)/dt=ω2
!d(ω1)/dt=[      ss1 - L21*μ2*cosΔ*ss2 ] /D
!d(ω2)/dt=[ -cosΔ/L21*ss1 +        ss2 ] /D
!
!---------- 微分方程式のまま、ルンゲ・クッタ法で描画。
LET θ1=PI*0.8 ! 初期角度
LET θ2=PI*0.8
LET w1=0 ! 初期角速度
LET w2=0
!
DEF ss1(w2,θ1,θ2)=-g/L1*SIN(θ1) -μ2*L21*w2^2*SIN(θ1-θ2)
DEF ss2(w1,θ1,θ2)=-g/L2*SIN(θ2) +        w1^2*SIN(θ1-θ2)/L21
DEF D(θ1,θ2)=1-μ2*COS(θ1-θ2)^2
!
DEF α1(w1,w2,θ1,θ2)=( ss1(w2,θ1,θ2) -L21*μ2*COS(θ1-θ2)*ss2(w1,θ1,θ2) )/D(θ1,θ2)
DEF α2(w1,w2,θ1,θ2)=(-ss1(w2,θ1,θ2)*COS(θ1-θ2)/L21     +ss2(w1,θ1,θ2) )/D(θ1,θ2)

SUB RungeKutta
   LET w11=w1
   LET w12=w2
   LET α11=α1(w1,w2,θ1,θ2)
   LET α12=α2(w1,w2,θ1,θ2)
   !
   LET w21=w1+α11*dt/2
   LET w22=w2+α12*dt/2
   LET α21=α1(w21,w22,θ1+w11*dt/2,θ2+w12*dt/2)
   LET α22=α2(w21,w22,θ1+w11*dt/2,θ2+w12*dt/2)
   !
   LET w31=w1+α21*dt/2
   LET w32=w2+α22*dt/2
   LET α31=α1(w31,w32,θ1+w21*dt/2,θ2+w22*dt/2)
   LET α32=α2(w31,w32,θ1+w21*dt/2,θ2+w22*dt/2)
   !
   LET w41=w1+α31*dt
   LET w42=w2+α32*dt
   LET α41=α1(w41,w42,θ1+w31*dt,θ2+w32*dt)
   LET α42=α2(w41,w42,θ1+w31*dt,θ2+w32*dt)
   !
   LET θ1=θ1+(w11+2*w21+2*w31+w41)*dt/6
   LET θ2=θ2+(w12+2*w22+2*w32+w42)*dt/6
   LET w1=w1+(α11+2*α21+2*α31+α41)*dt/6
   LET w2=w2+(α12+2*α22+2*α32+α42)*dt/6
END SUB

!----エネルギー・メーター
DEF ep1=m1*g*(L1-L1*COS(θ1)) !位置1
DEF em1=(L1*w1)^2*m1/2 !運動1
DEF ep2=m2*g*( (L1-L1*COS(θ1)+L2)-L2*COS(θ2) ) !位置2
DEF em2=( (L1*w1)^2+(L2*w2)^2-2*L1*w1*L2*w2*COS(PI+θ1-θ2) )*m2/2 !運動2
!
!----run
LET w=13
SET WINDOW -w,w,-w,w
SET COLOR MIX(15) .5,.5,.5
SET TEXT background "OPAQUE"
LET t0=TIME
DO
   LET t=TIME
   IF dt=<ABS(t-t0) THEN
      SET DRAW mode hidden
      CLEAR
      DRAW grid(5,5)
      PLOT TEXT,AT -w*0.92,w*0.9:"おもりのエネルギー[J]"
      PLOT TEXT,AT -w*0.92,w*0.83:"位置1 運動1  位置2 運動2"
      PLOT TEXT,AT -w*0.96,w*0.76,USING"##.#### ##.#### ##.#### ##.####":ep1,em1,ep2,em2
      PLOT TEXT,AT -w*0.86,w*0.69,USING"##.####     ##.####":ep1+em1,ep2+em2
      PLOT TEXT,AT -w*0.62,w*0.62,USING"##.####":ep1+em1+ep2+em2
      PLOT TEXT,AT w*0.25,w*0.9:"マウス 右ボタンで、終了。"
      PLOT TEXT,AT w*0.4,w*0.76,USING"演算ピッチ=#.### 秒":dt
      PLOT TEXT,AT w*0.4,w*0.69,USING"描画ピッチ=#.### 秒":t-t0
      LET t0=t
      DRAW Pendulum0 WITH ROTATE(θ1)
      CALL RungeKutta ! 次のθ1,θ2 へ更新
      SET DRAW mode explicit
      !stop
   END IF
   WAIT DELAY 0 ! 省電力効果
   MOUSE POLL mx,my,mlb,mrb
LOOP UNTIL mrb=1

PICTURE Pendulum0
   DRAW circle WITH SCALE(0.2)
   DRAW Pendulum1(L1,r1,"1")
   DRAW Pendulum1(L2,r2,"2") WITH ROTATE(θ2-θ1)*SHIFT(0,-L1)
END PICTURE

PICTURE Pendulum1(L,r,w$)
   PLOT LINES: 0,0;0,-L
   DRAW disk WITH SCALE(r)*SHIFT(0,-L)
   PLOT TEXT,AT r, r-L:w$
END PICTURE

END
 

戻る