本章首先介紹單擺和雙擺系統(tǒng)的公式推導(dǎo)，然后通過(guò)odeint()對(duì)其進(jìn)行數(shù)值求解并制作動(dòng)畫(huà)演示程序。

18.1  單擺模擬

如圖18-1所示，有一根不可伸長(zhǎng)、質(zhì)量不計(jì)的細(xì)棒，上端固定，下端系一質(zhì)點(diǎn)，這樣的裝置叫做單擺。

根據(jù)牛頓力學(xué)定律，我們可以列出如下微分方程：

其中，θ為單擺的擺角，l為單擺的長(zhǎng)度，g為重力加速度。

此微分方程的符號(hào)解無(wú)法直接求出，因此只能調(diào)用odeint()對(duì)其求數(shù)值解。

odeint()的調(diào)用參數(shù)如下：

odeint(func, y0, t, ...)

其中，func是Python的一個(gè)函數(shù)對(duì)象，用來(lái)計(jì)算微分方程組中每個(gè)未知函數(shù)的導(dǎo)數(shù)；y0為微分方程組中每個(gè)未知函數(shù)的初始值；t為需要進(jìn)行數(shù)值求解的時(shí)間點(diǎn)。它返回的是一個(gè)二維數(shù)組result，其第0軸的長(zhǎng)度為t的長(zhǎng)度，第1軸的長(zhǎng)度為變量的個(gè)數(shù)，因此result[:,i]為第i個(gè)未知函數(shù)的解。

計(jì)算微分的func函數(shù)的調(diào)用參數(shù)為func(y, t)，其中y是一個(gè)數(shù)組，為每個(gè)未知函數(shù)在t時(shí)刻的值，而func的返回值是每個(gè)未知函數(shù)在t時(shí)刻的導(dǎo)數(shù)。

odeint()要求每個(gè)微分方程只包含一階導(dǎo)數(shù)，因此我們需要對(duì)前面的微分方程進(jìn)行如下變形：