TI-Nspire版Pythonによる振り子のシミュレーション / 二重振り子

https://qiita.com/cyq04000/items/74541609ac03dbc598b9 の続き    

ここまでをまとめた _pend_simulator.tns_ は https://github.com/ti-nspire/tns_files/tree/master/pend_simulator に置いた。       

#概要
今度は二重振り子をシミュレートしてみる。

#サンプルファイル

```example_pend_double.py
from PEND_DOUBLE_Class import *
from CREATE_SIM_SPACE_Class import *
import ti_draw as td
import math

# 画面の大きさを決める。
WIDTH, HIGHT = td.get_screen_dim()
ASPECT = WIDTH / HIGHT
YMIN = -2.2
YMAX = -YMIN
XMIN = YMIN * ASPECT
XMAX = YMAX * ASPECT

# シミュレートすべきオブジェクトをいくつか作る。
"""
(内側振り子長さ,内側振り子角度,内側振り子速度,内側錘質量,内側錘密度, 
 外側振り子長さ,外側振り子角度,外側振り子速度,外側錘質量,外側錘密度,
1ステップ,1ステップ分割数)
"""
h = 1/2**5
num_of_div = 2**1
p0 = PEND_DOUBLE_Class(1, math.pi*0.999, 0, 2, 30, 1, math.pi*0.999, 0, 2, 30, h, num_of_div)
p1 = PEND_DOUBLE_Class(1.0001, math.pi*0.999, 0, 2, 30, 1, math.pi*0.999, 0, 2, 30, h, num_of_div)

# オブジェクトをシミュレーション空間に投入し、シミュレーションを始める。
CREATE_SIM_SPACE_Class([p0, p1], XMIN, XMAX, YMIN, YMAX).start_sim()
```

##実行結果
[![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/193936/c4d34635-4bfb-5efd-c897-45d8b6e71cd6.png)](https://youtu.be/sliPt0zidg8)

##二重振り子の質点のxy座標を計算するためのクラス

<details><summary>`ソースコード (PEND_DOUBLE_Class.py) を見る/隠す`</summary><div>

```PEND_DOUBLE_Class.py
from RK4_Class import *
import time
import ti_draw as td
import math

class PEND_DOUBLE_Class:
  def __init__(self, len0, rad0, vel0, m0, dens0, len1, rad1,vel1, m1, dens1, h=1/2**5, numOfDiv=1):

    vol0= m0 / dens0 # 重りの体積〔m^3〕。
    vol1= m1 / dens1 # 重りの体積〔m^3〕。
    self.MASS_R0 = ((3 * vol0) / (4 * math.pi))**(1 / 3) # 重りが球体の場合の半径〔m〕。
    self.MASS_R1 = ((3 * vol1) / (4 * math.pi))**(1 / 3) # 重りが球体の場合の半径〔m〕。

    self.LENGTH0 = len0
    self.LENGTH1 = len1

    # リセット用に引数を取っておく。
    self.ARGs = [len0, rad0, vel0, m0, dens0, len1, rad1,vel1, m1, dens1, h, numOfDiv]

    #更新すべき初期値を設定する。
    self.t0 = 0
    self.rad0 = rad0
    self.vel0 = vel0
    self.rad1 = rad1
    self.vel1 = vel1
    self.inits = [self.rad0, self.vel0, self.rad1, self.vel1]

    # 解くべき微分方程式を定義する。
    def rad0Dot(t, rad0, vel0, rad1, vel1):
      return vel0
    def vel0Dot(t, rad0, vel0, rad1, vel1):
      return -(9.80665*(m0*math.sin(rad0)-m1*(math.cos(rad0-rad1)*math.sin(rad1)-math.sin(rad0)))+(len1*vel1**2*math.cos(rad0-rad1)+len1*vel1**2)*m1*math.sin(rad0-rad1))/(len0*(m0+m1*(math.sin(rad0-rad1))**2)) 
    def rad1Dot(t, rad0, vel0, rad1, vel1):
      return vel1
    def vel1Dot(t, rad0, vel0, rad1, vel1):
      return (9.80665*(m0+m1)*(math.sin(rad0)*math.cos(rad0-rad1)-math.sin(rad1))+(len0*(m0+m1)*vel0**2+len1*m1*vel1**2*math.cos(rad0-rad1))*math.sin(rad0-rad1))/(len1*(m0+m1*(math.sin(rad0-rad1))**2))
    self.FUNCS = [rad0Dot, vel0Dot, rad1Dot, vel1Dot]

    # 上記の微分方程式を上記の初期値で解くための函数を実体化する。
    self.sol = RK4_Class(self.FUNCS, self.inits, h, numOfDiv)

  def update(self):
    self.sol.update() # 初期値を更新して、
    self.t0 = self.sol.t0 # 更新された時間を取り出して、
    self.rad0 = self.sol.inits[0] # 更新された角度(内側振り子)を取り出して、
    self.rad1 = self.sol.inits[2] # 更新された角度(外側振り子)を取り出して、

    return self


  def draw(self):
    x0 = self.LENGTH0 * math.sin(self.rad0) # 角度から内側質点のx座標を計算する。
    y0 = self.LENGTH0 * -math.cos(self.rad0) # 角度から内側質点のy座標を計算する。
    x1 = x0 + self.LENGTH1 * math.sin(self.rad1) # 角度から外側質点のx座標を計算する。
    y1 = y0 + self.LENGTH1 * -math.cos(self.rad1) # 角度から外側質点のy座標を計算する。
    
    td.draw_poly([0, x0, x1], [0, y0, y1]) # 原点からその質点まで線を引く。

    td.draw_circle(x0, y0, self.MASS_R0) # 質点の位置に円を描く。
    td.draw_circle(x1, y1, self.MASS_R1) # 質点の位置に円を描く。

    return self

  def _print(self):
    print(self.t0, self.sol.inits) # 計算結果の確認用。
    #print(self.t0, self.rad0, self.rad1) # 計算結果の確認用。

########
# test #
########
#if __name__ == "__main__":
def main():
  # len0, rad0, vel0, m0, dens0, len1, rad1,vel1, m1, dens1, h=1/2**5, numOfDiv=1
  pend = PEND_DOUBLE_Class(1, 3, 0, 2, 10, 1, 3, 0, 2, 10, 1/2**5, 2**0)
  for i in range(10):
    pend.update()._print()
```