# 雙擺系統的動畫模擬
相關文檔: [_單擺和雙擺模擬_](double_pendulum.html)
 
## 用odeint解雙擺系統
文件名: double_pendulum_odeint.py
```
# -*- coding: utf-8 -*-
from math import sin,cos
import numpy as np
from scipy.integrate import odeint
g = 9.8
class DoublePendulum(object):
def __init__(self, m1, m2, l1, l2):
self.m1, self.m2, self.l1, self.l2 = m1, m2, l1, l2
self.init_status = np.array([0.0,0.0,0.0,0.0])
def equations(self, w, t):
"""
微分方程公式
"""
m1, m2, l1, l2 = self.m1, self.m2, self.l1, self.l2
th1, th2, v1, v2 = w
dth1 = v1
dth2 = v2
#eq of th1
a = l1*l1*(m1+m2) # dv1 parameter
b = l1*m2*l2*cos(th1-th2) # dv2 paramter
c = l1*(m2*l2*sin(th1-th2)*dth2*dth2 + (m1+m2)*g*sin(th1))
#eq of th2
d = m2*l2*l1*cos(th1-th2) # dv1 parameter
e = m2*l2*l2 # dv2 parameter
f = m2*l2*(-l1*sin(th1-th2)*dth1*dth1 + g*sin(th2))
dv1, dv2 = np.linalg.solve([[a,b],[d,e]], [-c,-f])
return np.array([dth1, dth2, dv1, dv2])
def double_pendulum_odeint(pendulum, ts, te, tstep):
"""
對雙擺系統的微分方程組進行數值求解,返回兩個小球的X-Y坐標
"""
t = np.arange(ts, te, tstep)
track = odeint(pendulum.equations, pendulum.init_status, t)
th1_array, th2_array = track[:,0], track[:, 1]
l1, l2 = pendulum.l1, pendulum.l2
x1 = l1*np.sin(th1_array)
y1 = -l1*np.cos(th1_array)
x2 = x1 + l2*np.sin(th2_array)
y2 = y1 - l2*np.cos(th2_array)
pendulum.init_status = track[-1,:].copy() #將最后的狀態賦給pendulum
return [x1, y1, x2, y2]
if __name__ == "__main__":
import matplotlib.pyplot as pl
pendulum = DoublePendulum(1.0, 2.0, 1.0, 2.0)
th1, th2 = 1.0, 2.0
pendulum.init_status[:2] = th1, th2
x1, y1, x2, y2 = double_pendulum_odeint(pendulum, 0, 30, 0.02)
pl.plot(x1,y1, label = u"上球")
pl.plot(x2,y2, label = u"下球")
pl.title(u"雙擺系統的軌跡, 初始角度=%s,%s" % (th1, th2))
pl.legend()
pl.axis("equal")
pl.show()
```
## 擺動動畫
文件名: double_pendulum_animation.py
```
# -*- coding: utf-8 -*-
import matplotlib
matplotlib.use('WXAgg') # do this before importing pylab
import matplotlib.pyplot as pl
from double_pendulum_odeint import double_pendulum_odeint, DoublePendulum
fig = pl.figure(figsize=(4,4))
line1, = pl.plot([0,0], [0,0], "-o")
line2, = pl.plot([0,0], [0,0], "-o")
pl.axis("equal")
pl.xlim(-4,4)
pl.ylim(-4,2)
pendulum = DoublePendulum(1.0, 2.0, 1.0, 2.0)
pendulum.init_status[:] = 1.0, 2.0, 0, 0
x1, y1, x2, y2 = [],[],[],[]
idx = 0
def update_line(event):
global x1, x2, y1, y2, idx
if idx == len(x1):
x1, y1, x2, y2 = double_pendulum_odeint(pendulum, 0, 1, 0.05)
idx = 0
line1.set_xdata([0, x1[idx]])
line1.set_ydata([0, y1[idx]])
line2.set_xdata([x1[idx], x2[idx]])
line2.set_ydata([y1[idx], y2[idx]])
fig.canvas.draw()
idx += 1
import wx
id = wx.NewId()
actor = fig.canvas.manager.frame
timer = wx.Timer(actor, id=id)
timer.Start(1)
wx.EVT_TIMER(actor, id, update_line)
pl.show()
```
- 用Python做科學計算
- 軟件包的安裝和介紹
- NumPy-快速處理數據
- SciPy-數值計算庫
- matplotlib-繪制精美的圖表
- Traits-為Python添加類型定義
- TraitsUI-輕松制作用戶界面
- Chaco-交互式圖表
- TVTK-三維可視化數據
- Mayavi-更方便的可視化
- Visual-制作3D演示動畫
- OpenCV-圖像處理和計算機視覺
- Traits使用手冊
- 定義Traits
- Trait事件處理
- 設計自己的Trait編輯器
- Visual使用手冊
- 場景窗口
- 聲音的輸入輸出
- 數字信號系統
- FFT演示程序
- 頻域信號處理
- Ctypes和NumPy
- 自適應濾波器和NLMS模擬
- 單擺和雙擺模擬
- 分形與混沌
- 關于本書的編寫
- 最近更新
- 源程序集
- 三角波的FFT演示
- 在traitsUI中使用的matplotlib控件
- CSV文件數據圖形化工具
- NLMS算法的模擬測試
- 三維標量場觀察器
- 頻譜泄漏和hann窗
- FFT卷積的速度比較
- 二次均衡器設計
- 單擺擺動周期的計算
- 雙擺系統的動畫模擬
- 繪制Mandelbrot集合
- 迭代函數系統的分形
- 繪制L-System的分形圖