概要

TI-NspireでMicroPythonがサポートされた。ためしに単振り子の運動方程式(微分方程式)を近似的に解いて質点の動きをアニメ化してみる。

リアルタイムシミュレーションにしたいので、現状を描画 → 1ステップだけ未来を計算 → その状態を描画 → さらに1ステップだけ未来を計算 → ... を繰り返すことにする。

サンプルファイル

example_pend_single.py

from PEND_SINGLE_Class import *
from CREATE_SIM_SPACE_Class import *
from math import pi as PI

# シミュレートすべきオブジェクトをいくつか作る。
# 引数は(振り子長さ,錘描画半径,初期角度,初期角速度,1ステップ<,1ステップ分割数>)。
# 角度は、鉛直方向下向きを0ラジアン、反時計回りをプラス、時計回りをマイナスとする。
# 錘描画半径はシミュレーションには無関係。
h = 1/2**5
num_of_div = 2**2
p0 = PEND_SINGLE_Class(1, 0.1, PI * 0.998986328504169326 , 0, h, num_of_div)
p1 = PEND_SINGLE_Class(0.5, 0.1, PI * 0.7, 0, h, num_of_div)
p2 = PEND_SINGLE_Class(0.2, 0.1, PI * 0.5, 0, h, num_of_div)

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

実行結果

ハンドヘルドでもそこそこ動いている。

常微分方程式を解くためのクラス

Nspireにはrk23()という組み込み函数が用意されているが、ここでは古典的なルンゲクッタ法で解くことにする。外部モジュールは何も使っていないので普通のPythonでそのまま動く。

`ソースコード (RK4_Class.py) を見る/隠す`

RK4_Class.py

# Python による常微分方程式の数値解法_改良版 / 古典的 Runge-Kutta 法
class RK4_Class:
  def __init__(self, funcs, inits, h, num_of_div=1):
    _num_of_funcs = len(funcs)
    
    self.funcs = funcs
    self.t0 = 0
    self.inits = inits
    self.h = h / num_of_div
    self.half_h = self.h / 2
    self.sixth_h = self.h / 6
    self.range_dim = range(_num_of_funcs)
    self.range_div = range(num_of_div)
    self.f = [[0] * _num_of_funcs] * 4
    self.temp = [[0] * _num_of_funcs] * 3

  def update(self):
    for i in self.range_div:
      self.f[0] = [self.funcs[j](self.t0, *self.inits) for j in self.range_dim]
      self.temp[0] = [self.inits[j] + self.half_h * self.f[0][j] for j in self.range_dim]
      self.f[1] = [self.funcs[j](self.t0 + self.half_h, *self.temp[0]) for j in self.range_dim]
      self.temp[1] = [self.inits[j] + self.half_h * self.f[1][j] for j in self.range_dim]
      self.f[2] = [self.funcs[j](self.t0 + self.half_h, *self.temp[1]) for j in self.range_dim]
      self.temp[2] = [self.inits[j] + self.h * self.f[2][j] for j in self.range_dim]
      self.f[3] = [self.funcs[j](self.t0 + self.h, *self.temp[2]) for j in self.range_dim]
      self.inits = [self.inits[j] + self.sixth_h * (self.f[0][j] + 2 * (self.f[1][j] + self.f[2][j]) + self.f[3][j]) for j in self.range_dim]
      self.t0 += self.h
    return self

  def print(self):
    print(self.t0, self.inits)
    return self

########
# test #
########
# if __name__ == "__main__":
def main():
  # 解くべき連立常微分方程式を定義する。
  # 確認のため、厳密解の存在する方程式を定義する。単振り子とは無関係。
  def xDot(t, x, y): # x'(t) = y(t)
    return y
  def yDot(t, x, y): # y'(t) = t - x(t)
    return t - x
  funcs = [xDot, yDot]

  # 従属変数の初期値を指定する。
  x0 = 0
  y0 = 0
  inits = [x0, y0]

  # ステップ幅を指定する。
  h = 0.25 # =1/2**2
  
  # 1ステップだけ計算する函数を実体化して、
  rk4 = RK4_Class(funcs, inits, h)

  # 初期値を何度か更新して確認する。
  while rk4.t0 < 7.5:
    rk4.update().print()

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

`ソースコード (PEND_SINGLE_Class.py) を見る/隠す`

PEND_SINGLE_Class.py

from RK4_Class import *
import time
import ti_draw as td
import math

PI = math.pi

class PEND_SINGLE_Class:
  # (長さ,錘描画半径,初期角度,初期速度,1ステップ,1ステップ内部分割数)
  def __init__(self, length, mass_r, init_angle, init_vel, h=1/2**5, numOfDiv=1):
    
    # リセット用に引数を取っておく。
    self.ARGs = [length, mass_r, init_angle, init_vel, h, numOfDiv]

    self.LENGTH = length
    self.MASS_R = mass_r

    #更新すべき初期値を設定する。
    self.t0 = 0
    self.angle = init_angle
    self.vel = init_vel
    self.inits = [self.angle, self.vel]

    # 解くべき微分方程式を定義する。
    def dot_angle(t, angle, vel):
      return vel
    def dot_vel(t, angle, vel):
      return (-9.80665 / self.LENGTH) * math.sin(angle)
    self.FUNCS = [dot_angle, dot_vel]

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

  def update(self):
    self.rk4.update() # 初期値を更新して、
    self.t0 = self.rk4.t0 # 更新された時間を取り出して、
    self.angle = self.rk4.inits[0] # 更新された角度を取り出して、
    self.vel = self.rk4.inits[1] # 更新された角速度を取り出す。
    return self

  def draw(self):
    x = self.LENGTH * math.sin(self.angle) # 角度から質点のx座標を計算する。
    y = self.LENGTH * -math.cos(self.angle) # 角度から質点のy座標を計算する。
    td.draw_line(0, 0, x, y) # 原点からその質点まで線を引く。
    td.draw_circle(x, y, self.MASS_R) # 質点の位置に円を描く。
    return self

  def _print(self):
    print(self.t0, self.angle, self.vel) # 計算結果の確認用。

########
# test #
########
# if __name__ == "__main__":
def main():
   pend = PEND_SINGLE_Class(1, 0.1, PI * 0.999, 0, 1/2**4, 4)
   for i in range(10):
     pend.update()._print()

シミュレーション空間を作ってシミュレートするためのクラス

`ソースコード (CREATE_SIM_SPACE_Class.py) を見る/隠す`

CREATE_SIM_SPACE_Class.py

from PEND_SINGLE_Class import *
import ti_draw as td
import ti_system as tsys
from math import pi as PI

class CREATE_SIM_SPACE_Class:
  def __init__(self, objects):
    WIDTH, HIGHT = td.get_screen_dim()
    ASPECT = WIDTH / HIGHT

    self.objects = objects
    self.MAX_LENGTH = max([object.LENGTH for object in objects])
    self.Y_MIN = -self.MAX_LENGTH * 1.1
    self.Y_MAX = self.MAX_LENGTH * 1.1
    self.X_MIN = self.Y_MIN * ASPECT
    self.X_MAX = self.Y_MAX * ASPECT
    td.set_window(self.X_MIN, self.X_MAX, self.Y_MIN, self.Y_MAX)

    self.LINE_SPACE = (self.Y_MAX - self.Y_MIN) / 10 
    self.MESSAGES = ["To reset, press the esc key.",
                     "To pause or resume, press the enter key.",
                     "To exit, press the F12 or Home key." ,]

  def _wait_for_pressing_enter(self):
    while tsys.get_key() != "enter":
      pass

  def _update_objects(self):
    for object in self.objects:
      object.update()

  def _draw_objects(self):
    for object in self.objects:
      object.draw()

  def _reset_objects(self):
    td.clear()
    for object in self.objects:
      object.__init__(*object.ARGs)
      object.draw()

  def _display_time(self):
     td.draw_text(self.X_MIN, self.Y_MAX - self.LINE_SPACE, "sec: "+str(self.objects[0].t0))

  def _display_messages(self):
     td.draw_text(self.X_MIN, self.Y_MIN + self.LINE_SPACE * 2.5, self.MESSAGES[0])
     td.draw_text(self.X_MIN, self.Y_MIN + self.LINE_SPACE * 1.5, self.MESSAGES[1])
     td.draw_text(self.X_MIN, self.Y_MIN + self.LINE_SPACE * 0.5, self.MESSAGES[2])

  def _pause_or_resume_objects(self):
    while True:
      key = tsys.get_key()
      if key == "esc":
        self._reset_objects()
        self._display_time()
        self._display_messages()
        td.paint_buffer()
      if key == "enter":
        break

  def start_sim(self):
    td.use_buffer()
    self._draw_objects() # 初期状態を描いて、
    self._display_time()
    self._display_messages()
    td.paint_buffer()

    self._wait_for_pressing_enter() # enterキーが押されるまで時間をとめて、
    while True:
      td.clear()
      self._update_objects() # 1ステップだけ未来を計算して、
      self._draw_objects() # 描き直す、を繰り返す。
      self._display_time()
      self._display_messages()
      td.paint_buffer()

      key = tsys.get_key()
      if key == "enter":
        self._pause_or_resume_objects()
      if key == "esc":
        self._reset_objects()
      #time.sleep_ms(20)


########
# test #
########
# if __name__ == "__main__":
def main():
  # (長さ,錘描画半径,初期角度,初期速度,1ステップ<,1ステップ内部分割数>)
  p = PEND_SINGLE_Class(1, 0.1, PI * 0.999 , 0, 1/2**4 )
  sute = CREATE_SIM_SPACE_Class([p,])
  sute.start_sim()