お願いです

 投稿者：sukehiroﾒｰﾙ  投稿日：2010年 1月24日(日)07時17分32秒
  私も十進BASICの愛用者です。
苦労してァグランジェの方程式と、MuPaDを駆使して、2重振り子の運動方程式c1'',
c2''を導き出し、動作シミュレーションをすることができました。

3重振り子に挑戦したのですが、収拾がつかなくなりました。
どなたか、プログラム作成し掲示していただけると嬉しいです 		 

Re: お願いです

 投稿者：山中和義  投稿日：2010年 1月25日(月)10時08分5秒
  > No.987[元記事へ]

sukehiroさんへのお返事です。
!直列多重振り子
!参考サイト　http://www.aihara.co.jp/~taiji/pendula-equations/present-node5.html

LET N=3 !振り子の数

DIM m(N) !振り子のおもりの質量
DATA 0.2, 0.3, 0.2
MAT READ m

DIM l(N) !振り子の糸の長さ
DATA 3, 3, 4
MAT READ l

DIM th(N) !振り子のy軸に対する角度
LET th(1)=2*PI/3 !初期値
LET th(2)=-PI/2
LET th(3)=PI/6

!θ[d]（d=1,n）に関するラグランジュ運動方程式
!   m[d,n]*l[d]^2*d/dt^2{θ[d]}
! + Σ[i=1,d-1] m[d,n]*l[d]*l[i]*d/dt^2{θ[i]}*C[d,i]
! + Σ[i=d+1,n] m[i,n]*l[d]*l[i]*d/dt^2{θ[i]}*C[d,i]
! + Σ[i=1,d-1] m[d,n]*l[d]*l[i]*(d/dt{θ[i])^2*S[d,i]
! + Σ[i=d+1,n] m[i,n]*l[d]*l[i]*(d/dt{θ[i])^2*S[d,i]
! + m[d,n]*g*l[d]*SIN(θ[d])
! = 0

LET G=9.8 !重力加速度

FUNCTION sm(k,l) !m[k,l]=Σ[j=k,l] mj
   local s,j
   LET s=0
   FOR j=k TO l
      LET s=s+m(j)
   NEXT j
   LET sm=s
END FUNCTION
DEF C(i,j)=COS(th(i)-th(j)) !C[i,j]=COS(θi-θj)
DEF S(i,j)=SIN(th(i)-th(j)) !S[i,j]=SIN(θi-θj)

!２階微分方程式を連立１階微分方程式にする
!　ωd=d/dt{θd}とすると、R*d/dt^2{θd}=vよりd/dt{ωd}=[INV(R)*v]d

DIM R(N,N),v(N)
SUB FNC(t,x(), f()) !連立微分方程式 d/dt{Xi}=f(t,Xi)
   FOR d=1 TO N !行
      FOR i=1 TO N !列
         IF i>d THEN !右上部分
            LET R(d,i)=sm(i,N)*l(d)*l(i)*C(d,i)
         ELSEIF i=d THEN !対角線
            LET R(d,i)=sm(d,N)*l(d)^2
         ELSE !左下部分
            LET R(d,i)=sm(d,N)*l(d)*l(i)*C(d,i)
         END IF
      NEXT i
   NEXT d
   !!!MAT PRINT R;

   FOR d=1 TO N !行
      LET ss=0
      FOR i=1 TO d-1
         LET ss=ss-sm(d,N)*l(d)*l(i)*x(i)^2*S(d,i)
      NEXT i
      FOR i=d+1 TO N
         LET ss=ss-sm(i,N)*l(d)*l(i)*x(i)^2*S(d,i)
      NEXT i
      LET ss=ss-sm(d,N)*G*l(d)*SIN(th(d))

      LET v(d)=ss
   NEXT d
   !!!MAT PRINT v;

   MAT R=INV(R)
   MAT f=R*v
END SUB


SET WINDOW -10,10,-12,8 !描画領域

DIM w(N) !振り子のy軸に対する角速度
MAT w=ZER

LET h=0.01 !時間刻み幅 ⊿t
FOR t=0 TO 30 STEP h !※調整すること

   SET DRAW mode hidden !ちらつき防止開始
   CLEAR
   DRAW PendulumN(1) !親から順に描画する
   SET DRAW mode explicit !ちらつき防止終了

   !WAIT DELAY 0.2 !※必要に応じて

   !４次のルンゲ・クッタ法で連立微分方程式 d/dt{Xi}=f(t,Xi)を解く
   DIM k1(N),k2(N),k3(N),k4(N), f(N)

   DIM x(N),TT(N)
   MAT x=w
   CALL FNC(t,x, f)
   MAT k1=h*f

   MAT TT=(1/2)*k1
   MAT x=w+TT
   CALL FNC(t+h/2,x, f)
   MAT k2=h*f

   MAT TT=(1/2)*k2
   MAT x=w+TT
   CALL FNC(t+h/2,x, f)
   MAT k3=h*f

   MAT x=w+k3
   CALL FNC(t+h,x, f)
   MAT k4=h*f

   FOR i=1 TO N
      LET w(i)=w(i)+(k1(i)+2*(k2(i)+k3(i))+k4(i))/6 !x(t+h)=x(t)+(k1+2*k2+2*k3+k4)/6
   NEXT i

   MAT TT=h*w !ωd=d/dt{θd}より
   MAT th=th+TT
   !!!MAT PRINT th;

NEXT t


PICTURE PendulumN(i) !ｎ重の振り子を描く
   LET LL=l(i)
   LET AA=th(i)
   DRAW Pendulum(LL,m(i)) WITH ROTATE(AA) !親
   IF i<N THEN DRAW PendulumN(i+1) WITH SHIFT(LL*SIN(AA),-LL*COS(AA)) !階層関係（子へ）
END PICTURE

PICTURE Pendulum(L,r) !原点を基準に振り子を描く
   PLOT LINES: 0,0; 0,-L !糸
   DRAW disk WITH SCALE(r)*SHIFT(0,-L) !おもり
END PICTURE

END
 

戻る