はじめに
この記事は、以前にProcessingで書いた記事をPythonで書き直したものです。そのためカオスの説明や用いた運動方程式については、こちらを参照してください。ここでは、用いたライブラリやコードを載せます。
インストールするライブラリ
pygame ver1.9.5
普通にインストールしても動くとは思いますが、動かない場合はpip install pygame==1.9.5等のコマンドでインストールしてみてください。
コード
実行するコード
n_double_pendulum.py
# -*- coding:utf-8 -*-
import pygame
from math import pi
import sys
from dp_class import dp
#windowの大きさ
window_size = 800
#フレームレート
FPS = 60
# define colors
WHITE = (255, 255, 255)
BLACK = (0, 0, 0)
#振り子の本数
n = 100
#ずらす角度
dtheta = 0.000001;
#二重振り子オブジェクトを作成
furiko = [dp(pi *0.95 ,pi * 0.99 + i * dtheta ,window_size / 2,window_size / 2,window_size,FPS) for i in range(n)]
def main():
pygame.init() #初期化
screen = pygame.display.set_mode((window_size,window_size)) # 画面の設定
pygame.display.set_caption("二重振り子") #タイトルバーに表示する名前
# フレームを管理する時計をclockに格納
clock = pygame.time.Clock()
while(1):
#フレームレート設定
clock.tick(FPS)
#画面を黒で塗りつぶし
screen.fill(BLACK)
#振り子を描画
for i in range(n):
#振り子の角度を計算
furiko[i].calc()
pygame.draw.line(screen, WHITE, (window_size / 2, window_size / 2), (furiko[i].x1,furiko[i].y1))
pygame.draw.line(screen, WHITE, (furiko[i].x1, furiko[i].y1), (furiko[i].x2,furiko[i].y2))
pygame.display.update() #画面を更新
#イベント処理
for event in pygame.event.get():
if event.type == pygame.locals.QUIT: #閉じるボタンが押されたら終了
pygame.quit() #Pygameの終了(画面閉じられる)
sys.exit()
if __name__ =="__main__":
main()
二重振り子クラス
dp_class.py
from math import sin,cos;
class dp():
#コンストラクタ
def __init__(self ,th1_ , th2_ , x_ , y_, window, FPS):
self.__th1 = th1_
self.__th2 = th2_
self.__dth1 = 0.0
self.__dth2 = 0.0
#fpsから計算の刻み幅決定
self.__dt = 1.0 / FPS
self.__window = window
self.__x = x_
self.__y = y_
#第一振子長さ、第二振子長さ,重り質量1,2,重力加速度
self.__r1 = 2
self.__r2 = 1
self.__m1 = 0.5
self.__m2 = 0.5
self.__g = 9.81
#ルンゲクッタの係数
self.__kt1 = [0.0] * 4
self.__kt2 = [0.0] * 4
self.__kd1 = [0.0] * 4
self.__kd2 = [0.0] * 4
#第一角速度を計算
def __Fth1(self , t1_, t2_, d1_, d2_):
return d1_
#第二角速度を計算
def __Fth2(self , t1_, t2_, d1_, d2_):
return d2_
#第一角加速度を計算
def __Fdth1(self , t1_, t2_, d1_, d2_):
dd1_ = (self.__g * self.__m1 * sin(t1_) + self.__g * self.__m2 * sin(t1_ - 2 * t2_) / 2.0 +self.__g * self.__m2 *sin(t1_) / 2.0 + self.__m2 * self.__r1 * sin(2 *t1_- 2 *t2_) * pow(d1_,2) / 2.0 + self.__m2 * self.__r2 * sin(t1_ - t2_) * pow(d2_ ,2)) / (self.__r1 *(-self.__m1 + self.__m2 * pow(cos(t1_ - t2_) , 2) -self.__m2))
return dd1_
#第二角加速祖を計算
def __Fdth2(self , t1_, t2_, d1_, d2_):
dd2_ = ((self.__m1 + self.__m2) * (self.__g * sin(t2_) - self.__r1 * sin(t1_ - t2_) *pow(d1_,2) ) - (self.__g * self.__m1 * sin(t1_) + self.__g * self.__m2 *sin(t1_) + self.__m2 * self.__r2 * sin(t1_ - t2_) * pow(d2_ , 2)) * cos(t1_ - t2_)) / (self.__r2 * (-self.__m1 + self.__m2 * pow(cos(t1_ - t2_),2) - self.__m2))
return dd2_
#ルンゲクッタ法を用いてリンク角度を計算
def __rk(self):
#第1次近似を計算
self.__kt1[0] = self.__dt * self.__Fth1( self.__th1, self.__th2, self.__dth1, self.__dth2)
self.__kd1[0] = self.__dt * self.__Fdth1( self.__th1, self.__th2, self.__dth1, self.__dth2)
self.__kt2[0] = self.__dt * self.__Fth2( self.__th1, self.__th2, self.__dth1, self.__dth2)
self.__kd2[0] = self.__dt * self.__Fdth2( self.__th1, self.__th2, self.__dth1, self.__dth2)
#第2次近似以降を計算
i = 1
while i <= 3:
if i != 3:
self.__kt2[i] = self.__dt * self.__Fth2(self.__th1 + self.__kt1[i-1] / 2,self.__th2 + self.__kt2[i-1] /2,self.__dth1 + self.__kd1[i-1] /2,self.__dth2 +self.__kd2[i-1] /2)
self.__kt1[i] = self.__dt * self.__Fth1(self.__th1 + self.__kt1[i-1] / 2,self.__th2 + self.__kt2[i-1] /2,self.__dth1 + self.__kd1[i-1] /2,self.__dth2 +self.__kd2[i-1] /2)
self.__kd1[i] = self.__dt * self.__Fdth1(self.__th1 + self.__kt1[i-1] / 2,self.__th2 + self.__kt2[i-1] /2,self.__dth1 + self.__kd1[i-1] /2,self.__dth2 +self.__kd2[i-1] /2)
self.__kd2[i] = self.__dt * self.__Fdth2(self.__th1 + self.__kt1[i-1] / 2,self.__th2 + self.__kt2[i-1] /2,self.__dth1 + self.__kd1[i-1] /2,self.__dth2 +self.__kd2[i-1] /2)
else:
self.__kt1[i] = self.__dt * self.__Fth1(self.__th1 + self.__kt1[i-1],self.__th2 + self.__kt2[i-1] ,self.__dth1 + self.__kd1[i-1] ,self.__dth2 +self.__kd2[i-1] )
self.__kt2[i] = self.__dt * self.__Fth2(self.__th1 + self.__kt1[i-1],self.__th2 + self.__kt2[i-1] ,self.__dth1 + self.__kd1[i-1] ,self.__dth2 +self.__kd2[i-1] )
self.__kd1[i] = self.__dt * self.__Fdth1(self.__th1 + self.__kt1[i-1],self.__th2 + self.__kt2[i-1] ,self.__dth1 + self.__kd1[i-1] ,self.__dth2 +self.__kd2[i-1] )
self.__kd2[i] = self.__dt * self.__Fdth2(self.__th1 + self.__kt1[i-1],self.__th2 + self.__kt2[i-1] ,self.__dth1 + self.__kd1[i-1] ,self.__dth2 +self.__kd2[i-1] )
i += 1
#dt後の座標を更新
self.__th1 += (self.__kt1[0] + 2 * self.__kt1[1] + 2 * self.__kt1[2] + self.__kt1[3]) /6.0
self.__th2 += (self.__kt2[0] + 2 * self.__kt2[1] + 2 * self.__kt2[2] + self.__kt2[3]) /6.0
self.__dth1 += (self.__kd1[0] + 2 * self.__kd1[1] + 2 * self.__kd1[2] + self.__kd1[3]) /6.0
self.__dth2 += (self.__kd2[0] + 2 * self.__kd2[1] + 2 * self.__kd2[2] + self.__kd2[3]) /6.0
#描画関数
def calc(self):
self.__rk()
#リンク端の座標を計算
self.x1 = self.__x + self.__r1 * self.__window/8 * sin(self.__th1)
self.y1 = self.__y + self.__r1 * self.__window/8 * cos(self.__th1)
self.x2 = self.x1 + self.__r2 * self.__window/8 * sin(self.__th2)
self.y2 = self.y1 + self.__r2 * self.__window/8 * cos(self.__th2)
上の二つのコード同じファイル内に保存して、n_double_pendulum.pyを実行してください。
アニメーション
$10^{-5}$[rad]だけずらした、100個の二重振り子のアニメーション
初期には変化のなかった振り子が、突然ばらばらに動き始めるます。
このことから初期値に微小な変化があったときに、長期的な予測が不可能になるというカオス的振る舞いを確認できました。
おわり
Pythonのゲーム用ライブラリPygameが使いやすかったので、以前に書いたものを書き直してみました。Processingのタグだとあまり見られないのでPythonに...
プログラム結果を、文字だけでなくアニメーションなどで表示できると楽しくプログラムを学べると思います。