２重振り子のカオス

!２重振り子のカオス
!OPTION ARITHMETIC NATIVE
OPTION ARITHMETIC DECIMAL

LET g=0.2 !0.3!0.5!1!9.8　overflow
LET m1=2
LET m2=2
LET L1=20
LET L2=20

!------------ ２重振り子の方程式
! 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)=-g/L1*SIN(θ1) -μ2*L21*w2^2*SIN(θ1-θ2)
DEF ss2(w1)=-g/L2*SIN(θ2) + w1^2*SIN(θ1-θ2)/L21
DEF D=1-μ2*COS(θ1-θ2)^2

DEF α1(w1,w2)=( ss1(w2) -L21*μ2*COS(θ1-θ2)*ss2(w1) )/D
DEF α2(w1,w2)=(-ss1(w2)*COS(θ1-θ2)/L21 +ss2(w1) )/D

SUB 次次角度
LET g11=α1(w1,w2)
LET g12=α2(w1,w2)
LET g21=α1(w1+g11/2,w2+g12/2)
LET g22=α2(w1+g11/2,w2+g12/2)
LET g31=α1(w1+g21/2,w2+g22/2)
LET g32=α2(w1+g21/2,w2+g22/2)
LET g41=α1(w1+g31,w2+g32)
LET g42=α2(w1+g31,w2+g32)

LET w1=w1+(g11+2*g21+2*g31+g41)/6
LET w2=w2+(g12+2*g22+2*g32+g42)/6
LET θ1=θ1+w1
LET θ2=θ2+w2
END SUB

!----run
SET WINDOW -60,60,-60,60
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 0,55:"マウス右ボタンで、終了。"
DRAW 振子0 WITH ROTATE(θ1)*SHIFT(x,y)
CALL 次次角度
SET DRAW mode explicit
!stop
END IF
WAIT DELAY 0.03
MOUSE POLL mx,my,mlb,mrb !マウスの状態取得。
LOOP UNTIL mrb=1

PICTURE 振子0
DRAW circle WITH SCALE(1)
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

２重振り子のカオス	返事を書く
SECOND 2006/12/06 02:13:41
!２重振り子のカオス !OPTION ARITHMETIC NATIVE OPTION ARITHMETIC DECIMAL LET g=0.2 !0.3!0.5!1!9.8　overflow LET m1=2 LET m2=2 LET L1=20 LET L2=20 !------------ ２重振り子の方程式 ! 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)=-g/L1SIN(θ1) -μ2L21w2^2SIN(θ1-θ2) DEF ss2(w1)=-g/L2SIN(θ2) + w1^2SIN(θ1-θ2)/L21 DEF D=1-μ2COS(θ1-θ2)^2 DEF α1(w1,w2)=( ss1(w2) -L21μ2COS(θ1-θ2)ss2(w1) )/D DEF α2(w1,w2)=(-ss1(w2)COS(θ1-θ2)/L21 +ss2(w1) )/D SUB 次次角度 LET g11=α1(w1,w2) LET g12=α2(w1,w2) LET g21=α1(w1+g11/2,w2+g12/2) LET g22=α2(w1+g11/2,w2+g12/2) LET g31=α1(w1+g21/2,w2+g22/2) LET g32=α2(w1+g21/2,w2+g22/2) LET g41=α1(w1+g31,w2+g32) LET g42=α2(w1+g31,w2+g32) LET w1=w1+(g11+2g21+2g31+g41)/6 LET w2=w2+(g12+2g22+2g32+g42)/6 LET θ1=θ1+w1 LET θ2=θ2+w2 END SUB !----run SET WINDOW -60,60,-60,60 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 0,55:"マウス右ボタンで、終了。" DRAW 振子0 WITH ROTATE(θ1)SHIFT(x,y) CALL 次次角度 SET DRAW mode explicit !stop END IF WAIT DELAY 0.03 MOUSE POLL mx,my,mlb,mrb !マウスの状態取得。 LOOP UNTIL mrb=1 PICTURE 振子0 DRAW circle WITH SCALE(1) 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