概要
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()