|
|
!２重振子メーター付（解析したい方へ）再投稿
!２重振り子特有の下降衝撃時に、その急変する過程が、
!演算ピッチの間に入って脱落し、全エネルギーが変動する計算エラーが見られる。
!演算ピッチが小さいと、緩和するが、低速パソコンでは、描画ピッチが伴わない。
!-----
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くらいが望ましいが、描画ピッチ(画面に表示)が、ついて来れなくなったら戻す。
! ２つのピッチが、ズレていると物理的な速度ではなくなる。
!
!------------ ２重振り子の方程式
! 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 共に、重力方向０からの左回り角）
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)) !位置１
DEF em1=(L1*w1)^2*m1/2 !運動１
DEF ep2=m2*g*( (L1-L1*COS(θ1)+L2)-L2*COS(θ2) ) !位置２
DEF em2=( (L1*w1)^2+(L2*w2)^2-2*L1*w1*L2*w2*COS(PI+θ1-θ2) )*m2/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:"位置１ 運動１ 位置２ 運動２"
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
|
|