3バネ結合振り子.bas

 投稿者:村上元佑  投稿日:2012年 8月19日(日)12時22分24秒
  こんな、実験してみました

'  3バネ結合振り子.bas  MicroSoft BASIC互換モード
'
'  ラグランジェ法
'  数式処理システムMuPadを使っての、式の整理仮定
'    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                 '質点mの軌跡表示する場合は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
'------------------------------------

 

戻る