|
|
こんな、実験してみました
' 3ﾊﾞﾈ結合振り子.bas MicroSoft BASIC互換ﾓｰﾄﾞ
'
' ﾗｸﾞﾗﾝｼﾞｪ法
' 数式処理ｼｽﾃﾑＭｕＰａｄを使っての、式の整理仮定
' reset();
' /* 各ﾊﾞﾈの長さ */
' l1:=( (x1-x)^2+(y1-y)^2 )^1/2;
' l2:=( (x2-x)^2+(y2-y)^2 )^1/2;
' l3:=( (x3-x)^2+(y3-y)^2 )^1/2;
'
' T:=1/2*m*(xd1^2+yd1^2); /* 運動エネルギー */
' U:=m*g*y+1/2*k1*(lk1-l1)^2+1/2*k2*(lk2-l2)^2+1/2*k3*(lk3-l3)^2; /* 位置エネルギー */
' L:=T-U; /* ラグラジアン */
'
' Lxd1:=diff(L,xd1); /* x'による偏微分dL/x' */
' Lyd1:=diff(L,yd1); /* y'による偏微分dL/y' */
'
' Lx:=diff(L,x); /* xによる偏微分dL/x */
' Ly:=diff(L,y); /* yによる偏微分dL/y */
'
' xd2:=Lx/m; /* x'' */
' yd2:=Ly/m; /* y'' */
'
' 下記はMuPadで求めたxd2(x''),yd2(y'')を、ﾌﾟﾛｸﾞﾗﾑ記述したもの
' Lx_1=-k1*(x1-x)*( lk1-(x1-x)^2/2-(y1-y)^2/2 )
' Lx_2=-k2*(x2-x)*( lk2-(x2-x)^2/2-(y2-y)^2/2 )
' Lx_3=-k3*(x2-x)*( lk3-(x3-x)^2/2-(y3-y)^2/2 )
' Lx=lx_1+lx_2+lx_3
' xd2=Lx/m
'
' Ly_1=-k1*(y1-y)*( lk1-(x1-x)^2/2-(y1-y)^2/2 )-g*m
' Ly_2=-k2*(y2-y)*( lk2-(x2-x)^2/2-(y2-y)^2/2 )
' Ly_3=-k3*(y3-y)*( lk3-(x3-x)^2/2-(y3-y)^2/2 )
' Ly=Ly_1+Ly_2+Ly_3
' yd2=Ly/m
'
'
'--------------------------------------------------------
width=1260*2*2
height=640
SET BITMAP SIZE width, height
reft=0
right=width
bottom=0
top=height
dx=0
dy=0
'左端，右端，下端，上端
set window left+dx,right,bottom+dy,top
'描画エリアの背景色着色範囲設定
set area color 1 !'黒
plot area : left+dx,bottom+dy;left+dx,top;right,top;right,bottom+dy
DRAW AXES (100,100) '(right/k,top/k)
xbase=200
ybase=100
'--------------------------------------------------------
flag=1 '質点ｍの軌跡表示する場合は１
x1=0:y1=500 'ﾊﾞﾈl1の支持点座標
x2=500+200:y2=500 'ﾊﾞﾈl2の支持点座標
x3=250:y3=0 'ﾊﾞﾈl3の支持点座標
x=300:y=250 '3つのﾊﾞﾈに接続した質点の初期座標
g=9.8 '重力定数
m=6 '質点質量
l1=( (x1-x)^2+(y1-y)^2 )^1/2 'ﾊﾞﾈ1の初期長さ
l2=( (x2-x)^2+(y2-y)^2 )^1/2 'ﾊﾞﾈ2の初期長さ
l3=( (x3-x)^2+(y3-y)^2 )^1/2 'ﾊﾞﾈ3の初期長さ
lk1=l1/2 'ﾊﾞﾈ1の、無張力時の長さ
lk2=l2/3 'ﾊﾞﾈ2の、無張力時の長さ
lk3=l3/4 'ﾊﾞﾈ3の、無張力時の長さ
k1=0.6 'ﾊﾞﾈ1の、ﾊﾞﾈ定数
k2=0.8 'ﾊﾞﾈ2の、ﾊﾞﾈ定数
k3=0.9 'ﾊﾞﾈ3の、ﾊﾞﾈ定数
dt=0.0002 '時間刻幅
call draw '描画ﾙｰﾁﾝ
input aa$
'--------------------------------------------------------
do while zz=0
Lx_1=-k1*(x1-x)*( lk1-(x1-x)^2/2-(y1-y)^2/2 )
Lx_2=-k2*(x2-x)*( lk2-(x2-x)^2/2-(y2-y)^2/2 )
Lx_3=-k3*(x3-x)*( lk3-(x3-x)^2/2-(y3-y)^2/2 )
Lx=Lx_1+Lx_2+Lx_3
xd2=Lx/m ' x''
Ly_1=-k1*(y1-y)*( lk1-(x1-x)^2/2-(y1-y)^2/2 )-g*m
Ly_2=-k2*(y2-y)*( lk2-(x2-x)^2/2-(y2-y)^2/2 )
Ly_3=-k3*(y3-y)*( lk3-(x3-x)^2/2-(y3-y)^2/2 )
Ly=Ly_1+Ly_2+Ly_3
yd2=Ly/m ' y''
xd1=xd1+xd2*dt '質点mの、x成分速度
x=x+xd1*dt '質点mの、x成分位置
yd1=yd1+yd2*dt '質点mの、y成分速度
y=y+yd1*dt '質点mの、y成分位置
call draw '描画ﾙｰﾁﾝ
loop
'------------------------------------
sub draw
'描画ﾙｰﾁﾝ
local r
r=10
set draw mode hidden
'表示消去
line(xbase+x1,ybase+y1)-(xbase+xo,ybase+yo),1 'l1の表示
line(xbase+x2,ybase+y2)-(xbase+xo,ybase+yo),1 'l2の表示
line(xbase+x3,ybase+y3)-(xbase+xo,ybase+yo),1 'l3の表示
circle(xbase+xo,ybase+yo),r,1 ',,,,f
'表示更新
line(xbase-100,ybase)-(1000,ybase),6 '下の線
line(xbase-100,ybase+500)-(1000,ybase+500),6 '上の線
line(xbase+x1,ybase+y1)-(xbase+x,ybase+y),5 'l1の表示
line(xbase+x2,ybase+y2)-(xbase+x,ybase+y),5 'l2の表示
line(xbase+x3,ybase+y3)-(xbase+x,ybase+y),5 'l3の表示
circle(xbase+x,ybase+y),r,4 ',,,,f
if flag=1 then
circle(xbase+x,ybase+y),2,4 '質点mの、軌跡表示
end if
set draw mode explicit
xo=x:yo=y
end sub
'------------------------------------
|
|