!修正、２重振り子

!修正、２重振り子
!　ルンゲ・クッタ法　部分を、修正した。

LET g=9.8
LET m1=0.005
LET m2=0.005
LET L1=1
LET L2=1

!------------ ２重振り子の方程式
! 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

!---------- 上式と同一。
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)

LET dt=0.02

SUB 次次角度
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(w1+α11*dt/2,w2+α12*dt/2,θ1+w11*dt/2,θ2+w12*dt/2)
LET α22=α2(w1+α11*dt/2,w2+α12*dt/2,θ1+w11*dt/2,θ2+w12*dt/2)
LET w31=w1+α21*dt/2
LET w32=w2+α22*dt/2
LET α31=α1(w1+α21*dt/2,w2+α22*dt/2,θ1+w21*dt/2,θ2+w22*dt/2)
LET α32=α2(w1+α21*dt/2,w2+α22*dt/2,θ1+w21*dt/2,θ2+w22*dt/2)
LET w41=w1+α31*dt
LET w42=w2+α32*dt
LET α41=α1(w1+α31*dt,w2+α32*dt,θ1+w31*dt,θ2+w32*dt)
LET α42=α2(w1+α31*dt,w2+α32*dt,θ1+w31*dt,θ2+w32*dt)
LET ww1=(w11+2*w21+2*w31+w41)/6
LET ww2=(w12+2*w22+2*w32+w42)/6
LET aa1=(α11+2*α21+2*α31+α41)/6
LET aa2=(α12+2*α22+2*α32+α42)/6

LET θ1=θ1+ww1*dt
LET θ2=θ2+ww2*dt
LET w1= w1+aa1*dt
LET w2= w2+aa2*dt
END SUB

!----run
LET w=3
SET WINDOW -w,w,-w,w
LET r1=SQR(m1)
LET r2=SQR(m2)
LET tb=999
LET t0=TIME
DO
LET t=(TIME-t0)
IF t<>tb THEN
LET tb=t
SET DRAW mode hidden
CLEAR
PLOT TEXT,AT w*0.25,w*0.9:"マウス右ボタンで、終了。"
DRAW 振子0 WITH ROTATE(θ1)*SHIFT(x,y)
CALL 次次角度
SET DRAW mode explicit
!stop
END IF
WAIT DELAY 0.01
MOUSE POLL mx,my,mlb,mrb !マウスの状態取得。
LOOP UNTIL mrb=1

PICTURE 振子0
DRAW circle WITH SCALE(w/80)
DRAW 振子1(L1,r1)
DRAW 振子1(L2,r2) WITH ROTATE(θ2-θ1)*SHIFT(0,-L1)
END PICTURE

PICTURE 振子1(L,r)
PLOT LINES: 0,0;0,-L
DRAW disk WITH SCALE(r)*SHIFT(0,-L)
END PICTURE

END

!更新。エネルギー・メーターで検査する。
!２重振り子特有の下降衝撃時に、全エネルギーが変動する計算エラーが見られる。
!演算ピッチが小さいと、緩和するが、低速パソコンでは、描画ピッチが伴わない。

!----２重振り子のカオス
LET g= 9.8 !m/s^2
LET m1=0.1 !kg
LET m2=0.1 !kg
LET L1= 5 !m
LET L2= 5 !m

LET dt=0.05 !sec. 演算ピッチ。高速パソコンほど、小さく。(0.05は、500MHz)

!※描画ピッチが、ついて来れなくなったら戻す。
!　２つのピッチが、ズレていると物理的な速度ではなくなる。

!------------ ２重振り子の方程式
! 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

!---------- 上式と同一。
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 次次角度
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=15
SET WINDOW -w,w,-w,w
LET r1=SQR(m1)
LET r2=SQR(m2)
LET t0=TIME
DO
LET t=TIME
IF dt=<ABS(t-t0) THEN
SET DRAW mode hidden
CLEAR
PLOT TEXT,AT w*0.25,w*0.9:"マウス右ボタンで、終了。"
PLOT TEXT,AT -w*0.9,w*0.9:"エネルギー[J]"
PLOT TEXT,AT -w*0.9,w*0.83:"位置１　運動１　位置２　運動２"
PLOT TEXT,AT -w*0.95,w*0.76,USING"##.#### ##.####　##.#### ##.####":ep1,em1,ep2,em2
PLOT TEXT,AT -w*0.85,w*0.69,USING"##.####　　　　　##.####":ep1+em1,ep2+em2
PLOT TEXT,AT -w*0.61,w*0.62,USING"##.####":ep1+em1+ep2+em2
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 振子0 WITH ROTATE(θ1)
CALL 次次角度
SET DRAW mode explicit
!stop
END IF
WAIT DELAY 0 !クロック・アップを押える。
MOUSE POLL mx,my,mlb,mrb !マウスの状態取得。
LOOP UNTIL mrb=1

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

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

END

!修正、２重振り子	返事を書く
SECOND 2006/12/16 08:17:59
!修正、２重振り子 !　ルンゲ・クッタ法　部分を、修正した。 LET g=9.8 LET m1=0.005 LET m2=0.005 LET L1=1 LET L2=1 !------------ ２重振り子の方程式 ! d(dθ1)/dt^2= ! [ g{sinθ2cosΔ-μsinθ1}-{L2(θ2/dt)^2+L1(θ1/dt)^2cosΔ}sinΔ] !　　　　　　　　　　/[ L1{μ-cosΔ^2}] ! d(dθ2)/dt^2= ! [ gμ{sinθ1cosΔ-sinθ2}+{μL1(θ1/dt)^2+L2(θ2/dt)^2cosΔ}sinΔ] !　　　　　　　　　　/[ L2{μ-cosΔ^2}] ! g= 重力加速度　μ=(m1+m2)/m2　Δ=θ1-θ2 !---------- 上式と同一。 LET μ2=m2/(m1+m2) LET L21=L2/L1 !ss1=-g/L1sinθ1 -μ2L21ω2^2sinΔ !ss2=-g/L2sinθ2 +　　　　ω1^2sinΔ/L21 !D=1-μ2COSΔ^2 !d(θ1)/dt=ω1 !d(θ2)/dt=ω2 !d(ω1)/dt=[　　　　　　ss1 - L21μ2cosΔss2 ] /D !d(ω2)/dt=[ -cosΔ/L21ss1 + 　　　　　　　ss2 ] /D !---------- 式は解けないまま、ルンゲ・クッタ法で描画。 LET θ1=PI0.8 !初期角度 LET θ2=PI0.8 LET w1=0 !初期角速度 LET w2=0 DEF ss1(w2,θ1,θ2)=-g/L1SIN(θ1) -μ2L21w2^2SIN(θ1-θ2) DEF ss2(w1,θ1,θ2)=-g/L2SIN(θ2) + w1^2SIN(θ1-θ2)/L21 DEF D(θ1,θ2)=1-μ2COS(θ1-θ2)^2 DEF α1(w1,w2,θ1,θ2)=( ss1(w2,θ1,θ2) -L21μ2COS(θ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) LET dt=0.02 SUB 次次角度 LET w11=w1 LET w12=w2 LET α11=α1(w1,w2,θ1,θ2) LET α12=α2(w1,w2,θ1,θ2) LET w21=w1+α11dt/2 LET w22=w2+α12dt/2 LET α21=α1(w1+α11dt/2,w2+α12dt/2,θ1+w11dt/2,θ2+w12dt/2) LET α22=α2(w1+α11dt/2,w2+α12dt/2,θ1+w11dt/2,θ2+w12dt/2) LET w31=w1+α21dt/2 LET w32=w2+α22dt/2 LET α31=α1(w1+α21dt/2,w2+α22dt/2,θ1+w21dt/2,θ2+w22dt/2) LET α32=α2(w1+α21dt/2,w2+α22dt/2,θ1+w21dt/2,θ2+w22dt/2) LET w41=w1+α31dt LET w42=w2+α32dt LET α41=α1(w1+α31dt,w2+α32dt,θ1+w31dt,θ2+w32dt) LET α42=α2(w1+α31dt,w2+α32dt,θ1+w31dt,θ2+w32dt) LET ww1=(w11+2w21+2w31+w41)/6 LET ww2=(w12+2w22+2w32+w42)/6 LET aa1=(α11+2α21+2α31+α41)/6 LET aa2=(α12+2α22+2α32+α42)/6 LET θ1=θ1+ww1dt LET θ2=θ2+ww2dt LET w1= w1+aa1dt LET w2= w2+aa2dt END SUB !----run LET w=3 SET WINDOW -w,w,-w,w LET r1=SQR(m1) LET r2=SQR(m2) LET tb=999 LET t0=TIME DO LET t=(TIME-t0) IF t<>tb THEN LET tb=t SET DRAW mode hidden CLEAR PLOT TEXT,AT w0.25,w0.9:"マウス右ボタンで、終了。" DRAW 振子0 WITH ROTATE(θ1)SHIFT(x,y) CALL 次次角度 SET DRAW mode explicit !stop END IF WAIT DELAY 0.01 MOUSE POLL mx,my,mlb,mrb !マウスの状態取得。 LOOP UNTIL mrb=1 PICTURE 振子0 DRAW circle WITH SCALE(w/80) DRAW 振子1(L1,r1) DRAW 振子1(L2,r2) WITH ROTATE(θ2-θ1)SHIFT(0,-L1) END PICTURE PICTURE 振子1(L,r) PLOT LINES: 0,0;0,-L DRAW disk WITH SCALE(r)SHIFT(0,-L) END PICTURE END

Re: !修正、２重振り子	返事を書く
SECOND 2006/12/17 01:23:11
!更新。エネルギー・メーターで検査する。 !２重振り子特有の下降衝撃時に、全エネルギーが変動する計算エラーが見られる。 !演算ピッチが小さいと、緩和するが、低速パソコンでは、描画ピッチが伴わない。 !----２重振り子のカオス LET g= 9.8 !m/s^2 LET m1=0.1 !kg LET m2=0.1 !kg LET L1= 5 !m LET L2= 5 !m LET dt=0.05 !sec. 演算ピッチ。高速パソコンほど、小さく。(0.05は、500MHz) !※描画ピッチが、ついて来れなくなったら戻す。 !　２つのピッチが、ズレていると物理的な速度ではなくなる。 !------------ ２重振り子の方程式 ! d(dθ1)/dt^2= ! [ g{sinθ2cosΔ-μsinθ1}-{L2(θ2/dt)^2+L1(θ1/dt)^2cosΔ}sinΔ] !　　　　　　　　　　/[ L1{μ-cosΔ^2}] ! d(dθ2)/dt^2= ! [ gμ{sinθ1cosΔ-sinθ2}+{μL1(θ1/dt)^2+L2(θ2/dt)^2cosΔ}sinΔ] !　　　　　　　　　　/[ L2{μ-cosΔ^2}] ! g= 重力加速度　μ=(m1+m2)/m2　Δ=θ1-θ2 !---------- 上式と同一。 LET μ2=m2/(m1+m2) LET L21=L2/L1 !ss1=-g/L1sinθ1 -μ2L21ω2^2sinΔ !ss2=-g/L2sinθ2 +　　　　ω1^2sinΔ/L21 !D=1-μ2COSΔ^2 !d(θ1)/dt=ω1 !d(θ2)/dt=ω2 !d(ω1)/dt=[　　　　　　ss1 - L21μ2cosΔss2 ] /D !d(ω2)/dt=[ -cosΔ/L21ss1 + 　　　　　　　ss2 ] /D !---------- 式は解けないまま、ルンゲ・クッタ法で描画。 LET θ1=PI0.8 !初期角度 LET θ2=PI0.8 LET w1=0 !初期角速度 LET w2=0 DEF ss1(w2,θ1,θ2)=-g/L1SIN(θ1) -μ2L21w2^2SIN(θ1-θ2) DEF ss2(w1,θ1,θ2)=-g/L2SIN(θ2) + w1^2SIN(θ1-θ2)/L21 DEF D(θ1,θ2)=1-μ2COS(θ1-θ2)^2 DEF α1(w1,w2,θ1,θ2)=( ss1(w2,θ1,θ2) -L21μ2COS(θ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 次次角度 LET w11=w1 LET w12=w2 LET α11=α1(w1,w2,θ1,θ2) LET α12=α2(w1,w2,θ1,θ2) LET w21=w1+α11dt/2 LET w22=w2+α12dt/2 LET α21=α1(w21,w22,θ1+w11dt/2,θ2+w12dt/2) LET α22=α2(w21,w22,θ1+w11dt/2,θ2+w12dt/2) LET w31=w1+α21dt/2 LET w32=w2+α22dt/2 LET α31=α1(w31,w32,θ1+w21dt/2,θ2+w22dt/2) LET α32=α2(w31,w32,θ1+w21dt/2,θ2+w22dt/2) LET w41=w1+α31dt LET w42=w2+α32dt LET α41=α1(w41,w42,θ1+w31dt,θ2+w32dt) LET α42=α2(w41,w42,θ1+w31dt,θ2+w32dt) LET θ1=θ1+(w11+2w21+2w31+w41)dt/6 LET θ2=θ2+(w12+2w22+2w32+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 !続く

Re: !更新。エネルギー・メーターで検査する。	返事を書く
SECOND 2006/12/17 01:24:39
!続く !----エネルギー・メーター DEF ep1=m1g(L1-L1COS(θ1)) !位置１ DEF em1=(L1w1)^2m1/2 !運動１ DEF ep2=m2g( (L1-L1COS(θ1)+L2)-L2COS(θ2) ) !位置２ DEF em2=( (L1w1)^2+(L2w2)^2-2L1w1L2w2COS(PI+θ1-θ2) )m2/2 !運動２ !----run LET w=15 SET WINDOW -w,w,-w,w LET r1=SQR(m1) LET r2=SQR(m2) LET t0=TIME DO LET t=TIME IF dt=<ABS(t-t0) THEN SET DRAW mode hidden CLEAR PLOT TEXT,AT w0.25,w0.9:"マウス右ボタンで、終了。" PLOT TEXT,AT -w0.9,w0.9:"エネルギー[J]" PLOT TEXT,AT -w0.9,w0.83:"位置１　運動１　位置２　運動２" PLOT TEXT,AT -w0.95,w0.76,USING"##.#### ##.####　##.#### ##.####":ep1,em1,ep2,em2 PLOT TEXT,AT -w0.85,w0.69,USING"##.####　　　　　##.####":ep1+em1,ep2+em2 PLOT TEXT,AT -w0.61,w0.62,USING"##.####":ep1+em1+ep2+em2 PLOT TEXT,AT w0.4,w0.76,USING"演算ピッチ=#.### 秒":dt PLOT TEXT,AT w0.4,w0.69,USING"描画ピッチ=#.### 秒":t-t0 LET t0=t DRAW 振子0 WITH ROTATE(θ1) CALL 次次角度 SET DRAW mode explicit !stop END IF WAIT DELAY 0 !クロック・アップを押える。 MOUSE POLL mx,my,mlb,mrb !マウスの状態取得。 LOOP UNTIL mrb=1 PICTURE 振子0 DRAW circle WITH SCALE(0.2) DRAW 振子1(L1,r1,"1") DRAW 振子1(L2,r2,"2") WITH ROTATE(θ2-θ1)SHIFT(0,-L1) END PICTURE PICTURE 振子1(L,r,w$) PLOT LINES: 0,0;0,-L DRAW disk WITH SCALE(r)*SHIFT(0,-L) PLOT TEXT,AT 0, 0.5-L:w$ END PICTURE END