20 拋物體

到目前為止，我們使用的微分方程是一階的，這意味著它們只涉及一階導(dǎo)數(shù)引几。在本章中，我們將注意力轉(zhuǎn)向二階常微分方程挽铁，它可能涉及一階和二階導(dǎo)數(shù)伟桅。

我們將重新訪問第1章中的“自由落體的硬幣”示例，并使用run_ode_solver分別在考慮和不考慮空氣阻力的情況下查找硬幣下降時的位置和速度屿储。

20.1 牛頓第二運動定律

一階常微分方程組可以寫作

其中 G 是 x 和 y 的某個函數(shù)（參見http://modsimpy.com/ode). 二階常微分方程組可寫作

其中 H 是 x , y 和 dy/dx 的函數(shù)贿讹。

在這一章里，我們將用到最著名并且很實用的牛頓第二運動定律够掠，其公式是一個二階常微分方程：

其中F是力或一組力的總和民褂，m 是運動物體的質(zhì)量，a 是它的加速度疯潭。

牛頓定律可能看起來不像微分方程赊堪，除非我們意識到加速度 a 是位置 y 對時間 t 的二階導(dǎo)數(shù)

所以牛頓第二定律也可以寫成

這絕對是二階常微分方程。一般來說竖哩，F 可以是時間哭廉、位置和速度的函數(shù)。

當(dāng)然相叁，這個“定律”實際上是一個模型遵绰，因為它是對現(xiàn)實世界的簡化。雖然這通常是正確的：

• 它只適用于 m 為常數(shù)的情況增淹。如果質(zhì)量取決于時間椿访、位置或速度，我們必須使用牛頓定律的更一般形式（見http://modsimpy.com/varmass).
• 對于非常微觀的事物虑润，它不是一個好的模型成玫，而另一個模型，量子力學(xué)可以更好地描述它們。
• 對于高速運動的物體哭当，它并不是一個很好的模型猪腕，可以用另一個模型，相對論力學(xué)來描述钦勘。

然而陋葡，對于質(zhì)量恒定、以中等速度運動的中型物體彻采，牛頓模型是非常有用的脖岛。如果我們能量化作用在這樣一個物體上的力，我們就可以預(yù)測它將如何移動颊亮。

20.2 拋硬幣

作為第一個例子，讓我們回想在1.1節(jié)中討論過的帝國大廈掉下來的硬幣陨溅。我們將實施兩種模式的系統(tǒng)：一種是不考慮空氣阻力终惑，另一種考慮空氣阻力。

考慮到帝國大廈高381米门扇，假設(shè)硬幣是從靜止?fàn)顟B(tài)掉落的雹有，初始條件是：

init = State(y=381 * m,
            v=0 * m/s)

其中 y 是地平面以上的高度，v 是速度臼寄。m 和 s 是 Pint 里提供的單位霸奕。

m = UNITS.meter
s = UNITS.second

唯一的系統(tǒng)參數(shù)是重力加速度：

g = 9.8 * m/s**2

此外，我們還將指定模擬的持續(xù)時間和步長：

t_end = 10 * s
dt = 0.1 * s

有了這些參數(shù)吉拳，時間步數(shù)是100质帅，這對于許多問題來說是很好的。一旦我們有了解決方案留攒，我們將增加步驟的數(shù)量煤惩，看看它對結(jié)果有什么影響。
我們需要一個系統(tǒng)來存儲參數(shù)：

system = System(init=init, g=g, t_end=t_end, dt=dt)

現(xiàn)在我們需要一個斜率函數(shù)炼邀，這就是問題的癥結(jié)所在魄揉。如我們所見，run_ode_solver可以解一階常微分方程組拭宁，但牛頓定律是二階常微分方程組洛退。但是，如果我們認識到這一點

1. 速度 v 是位置的導(dǎo)數(shù)杰标，dy/dt
2. 加速度 a 是速度的導(dǎo)數(shù)兵怯，dv/dt

我們可以將牛頓定律改寫為一階常微分方程組：

并且我們可以把這些等式轉(zhuǎn)換成一個斜率函數(shù)：

def slope_func(state, t, system):
    y, v = state
    
    dydt = v
    dvdt = -system.g
    
    return dydt, dvdt

第一個參數(shù)state包含硬幣的位置和速度。最后一個參數(shù)system包含系統(tǒng)參數(shù)g在旱，它是由于重力引起的加速度的大小摇零。

第二個參數(shù)t是時間。該斜率函數(shù)中不使用它，因為模型中的任何因素都不是時間相關(guān)的（見第18.1節(jié)）驻仅。我之所以包含它谅畅，是因為這個函數(shù)將由run_ode_solver調(diào)用，它總是提供相同的參數(shù)噪服，不管它們是否需要毡泻。

函數(shù)的其余部分是微分方程的直接轉(zhuǎn)換，替換a=-g粘优，這表示重力引起的加速度朝著y減小的方向發(fā)展仇味。slope_func返回一個包含兩個導(dǎo)數(shù)的序列。
在調(diào)用run_ode_solver之前雹顺，最好使用初始條件測試斜率函數(shù)：

dydt, dvdt = slpoe_func(system.init, 0, system)

20_1.png

圖20.1: 硬幣的高度與時間的關(guān)系丹墨，沒有空氣阻力

速度為 0 m/s，加速度為 9.8 m/s^2℃依ⅲ現(xiàn)在我們將run_ode_solver稱為如下：

results, details = run_ode_solver(system, slope_func)

結(jié)果是一個有兩列的時間框架(TimeFrame)：y 代表硬幣的高度贩挣；v 代表它的速度。

我們可以像這樣畫出結(jié)果圖：

def plot_position(results):
    plot(results.y)
    decorate(xlabel='Time (s)',
            ylabel='Position (m)')

圖20.1顯示了結(jié)果没酣。由于加速度恒定王财，速度線性增加，位置二次減少裕便；因此绒净，高度曲線是拋物線。

results.y的最后一個值是-109m偿衰，這意味著我們運行模擬的時間太長了挂疆。解決這個問題的一種方法是利用結(jié)果來估計硬幣落在地面上的時間。

ModSim庫提供交叉哎垦，它接受一個TimeSeries和一個值囱嫩，并在序列通過該值時返回一個時間序列。我們可以這樣找到硬幣高度為 0 的時間：

t_crossings = crossings(results.y, 0)

結(jié)果是一個只有一個值的數(shù)組漏设，8.818s∧校現(xiàn)在，我們可以用t_end=8.818再次運行模擬郑口，但是還有更好的方法鸳碧。

20.3 事件

作為一個選項，run_ode_solver可以使用一個event_function犬性，該函數(shù)檢測到一個“事件”瞻离，如硬幣擊中地面，那就結(jié)束模擬乒裆。

事件函數(shù)采用與斜率函數(shù)套利、state、t和system相同的參數(shù)。它們應(yīng)該返回一個在事件發(fā)生時變?yōu)?code>0的值肉迫。以下是一個事件函數(shù)验辞，用于檢測硬幣撞到地面上：

def event_func(state, t, system):
y, v = state
return y

返回值是硬幣的高度 y，當(dāng)硬幣碰到地面時喊衫，它變?yōu)?code>0跌造。

我們把事件函數(shù)像這樣傳遞給run_ode_solver：

results, details = run_ode_solver(system, slope_func,
                                  events=event_func)

然后我們可以得到在空中的時間和最終速度，如下所示：

t_sidewalk = get_last_label(results) * s
v_sidewalk = get_last_value(results.v)

在你繼續(xù)之前族购，你可能想讀一下本章的筆記本壳贪，chap20.ipynb，并練習(xí)一下寝杖。有關(guān)下載和運行代碼的說明违施，請參閱Section 0.4。

本書的中文翻譯由南開大學(xué)醫(yī)學(xué)院智能醫(yī)學(xué)工程專業(yè)2018級瑟幕、2019級的師生完成醉拓，方便后續(xù)學(xué)生學(xué)習(xí)《Python仿真建模》課程收苏。翻譯人員（排名不分前后）：薛淏源、金鈺愤兵、張雯鹿霸、張瑩睿、趙子雨秆乳、李翀懦鼠、慕振墺、許靖云屹堰、李文碩肛冶、尹瀛寰、沈紀辰扯键、迪力木拉睦袖、樊旭波、商嘉文荣刑、趙旭馅笙、連煦、楊永新厉亏、樊一諾董习、劉志鑫、彭子豪爱只、馬碧婷皿淋、吳曉玲、常智星、陳俊帆窝趣、高勝寒疯暑、韓志恒、劉天翔高帖、張藝瀟缰儿、劉暢。

整理校訂由劉暢完成散址，如果您發(fā)現(xiàn)有翻譯不當(dāng)或者錯誤乖阵，請郵件聯(lián)系changliu@nankai.edu.cn。

本書中文版不用于商業(yè)用途预麸，供大家自由使用瞪浸。

未經(jīng)允許，請勿轉(zhuǎn)載吏祸。