26_案例分析

26 案例學(xué)習(xí)

26.1 計算工具

在第20章，我們將一個二階微分方程改寫成一個一階方程，并且像這樣使用一個斜率函數(shù)解決方程：

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

我們使用crossings函數(shù)去尋找在仿真結(jié)果中的過零檢測。

然后我們用一個像這樣的事件函數(shù)：

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

當(dāng)一個事件發(fā)生時，就停止仿真。這時注意到，事件函數(shù)和斜率函數(shù)使用相同的參數(shù)。

在第21章，我們建立了一個空氣阻力模型并且使用了一個Params對象，這是一個參數(shù)集合：

params = Params(height = 381 * m,
                v_init = 0 * m / s,
                g = 9.8 * m/s**2,
                mass = 2.5e-3 * kg,
                diameter = 19e-3 * m,
                rho = 1.2 * kg/m**3,
                v_term = 18 * m / s)

那么，我們可以用一個新的方式創(chuàng)建一個System對象，就是復(fù)制來自Params對象的所有變量并進(jìn)行增加或改變變量：

    return System(params, area=area, C_d=C_d,
              init=init, t_end=t_end)

我們也可以用gradient函數(shù)去估計加速度和給定的速度：

a = gradient(results.v)

第22章引入了Vector對象，它們能夠在二維或三維空間中代表一些矢量，例如位移，速度，力和加速度。

A = Vector(3,4) * m

這一章還引入了軌跡圖，它可以在二維空間中顯示一個對象的路徑：

x = results.R.extract('x')
y = results.R.extract('y')

plot(x, y, label='trajectory')

在第23章，我們定義了一個可以計算棒球飛行距離的距離函數(shù)，它又可以作為發(fā)射角度的函數(shù)：

def range_func(angle, params):
    params = Params(params, angle=angle)
    system = make_system(params)
    results, details = run_ode_solver(system, slope_func,
                                       events=event_func)
    x_dist = get_last_value(results.R).x
    return x_dist

然后，我們用maximize尋找使距離最大的發(fā)射角度：

bounds = [0, 90] * degree
res = maximize(range_func, bounds, params)

這樣，你的工具包就完成了。第24章和第25章引入了旋轉(zhuǎn)物理學(xué)，但沒有新的計算工具。

26.2 蹦極

假設(shè)你想要創(chuàng)造“蹦極扣籃”的最高世界記錄，這是一個挑戰(zhàn)蹦極者在跳躍的最低點在一杯茶中泡曲奇的特技挑戰(zhàn)。在這個視頻里有一個例子：http://modsimpy.com/dunk.

既然紀(jì)錄是70米，那么我們設(shè)一個80米的蹦極目標(biāo)。我們會按照下面的模型假設(shè)開始：

  • 開始時，將蹦極繩索懸掛在一臺起重機(jī)上，連接點在一杯茶上方的80米處。

  • 在繩子完全伸展之前，它不會對蹦床者施加任何力量。這可能不是一個很好的假設(shè)；我們將重新開始。

  • 在繩子完全伸展之后，它遵循胡克定律；也就是說，它對線施加一個力，這個力與繩子超出其靜止長度的延伸成正比?？梢?http://modsimpy.com/hooke.

  • 蹦極者的質(zhì)量為75公斤。

  • 跳躍者受到阻力，所以他們的末速度為60米/秒。

我們的目標(biāo)是設(shè)置繩子的長度L，和它的彈簧系數(shù)k，這樣，跳躍著就會一直落到茶杯的位置并且不超過茶杯！

在這本書的資料里，你會找到一個筆記bungee.ipynb，它包含著本次案例學(xué)習(xí)的初始代碼和練習(xí)。

26.3 重寫蹦極扣籃

在前面的案例學(xué)習(xí)中，我們假設(shè)在繩被拉伸之前不會對跳躍者施加任何力。這是我們常說的因為繩子與跳躍者一起下落，所以繩子沒有作用，但這種直覺是不正確的。當(dāng)繩子下落時，它會將能量傳遞給跳躍者。

在http://modsimply.com/bungee，你會找到一篇文章[^1]解釋了這一過程，并得出跳躍者的加速度a，作為位移y和速度v的函數(shù)：

a = g + \frac{\mu v^2/2}{\mu(L+y)+2L}

g為重力加速度，L為繩子長度，\mu 為繩子質(zhì)量m和跳躍者質(zhì)量M的比。

如果你不相信模型的正確性，這個視頻會為你證實：http://modsimpy.com/drop.

在這本書的資料里，你會找到一個筆記bungee2.ipynb，它包含著本次案例學(xué)習(xí)的初始代碼和練習(xí)。當(dāng)我們改變繩子的質(zhì)量時，系統(tǒng)的行為會發(fā)生怎樣的變化？當(dāng)繩子的質(zhì)量等于跳躍者的質(zhì)量時，跳躍最低點的凈影響是什么？

[^1]Heck, Uylings, and Kdzierska, “Understanding the physics of bungee jumping”, Physics Education, Volume 45, Number 1, 2010.

26.4 蜘蛛俠

在這個案例研究中，我們將開發(fā)一個蜘蛛俠從連接在帝國大廈頂部的有彈性的織帶索上搖擺的模型。最開始，蜘蛛俠在鄰近的大廈頂部，就像圖26.1顯示的那樣。

26_1.png

圖26.1：蜘蛛俠案例研究的初始狀態(tài)圖。

原點O，在帝國大廈的底部。向量H代表織帶相對于O處連接大廈的位置。向量P代表蜘蛛俠相對于O處的位置。L表示從連接點到蜘蛛俠的向量。

根據(jù)下面的箭頭，從O到H，再到L，我們可以知道

H+L=P

那么我們可以像這樣計算L:

L=P-H

本次的案例學(xué)習(xí)目標(biāo)是：

1.實現(xiàn)預(yù)測蜘蛛俠軌跡的模型。

2.選擇合適的時間讓蜘蛛俠松開蛛絲，這樣他在著陸前能走得更遠(yuǎn)。

3.為蜘蛛俠選擇最好的角度跳下大樓，放開蛛網(wǎng)，最大限度地擴(kuò)大射程。

我們將使用以下參數(shù)：

1.根據(jù)蜘蛛俠的維基，蜘蛛俠重76千克。

2.假設(shè)他的末速度為60米/秒。

3.網(wǎng)的長度為100米。

4.網(wǎng)的初始角度是45度向左下方向。

5.當(dāng)繩子被拉伸時，蛛網(wǎng)的彈性系數(shù)為40N/m，當(dāng)它被壓縮時為0。

在這本書的資料里，你會找到一個筆記spiderman.ipynb，它包含著本次案例學(xué)習(xí)的初始代碼。通讀這個筆記并運(yùn)行代碼。代碼中使用了minimize，它是一個可以搜索最優(yōu)參數(shù)集的SciPy函數(shù)（與minimize_scalar相反，它只能沿著單個軸搜索）。

26.5 小貓

讓我們模擬一只小貓在鋪開廁紙。參考資料，請看這個視頻：http://modsimpy.com/kitten.

小貓和紙輥的相互作用是復(fù)雜的。為了簡單起見，我們假設(shè)小貓用恒定的力向下拉動紙輥的自由端。與此同時，我們將忽略滾輪與滾軸之間的摩擦。

圖26.2展示了帶有r，F(xiàn)和τ 的紙輥。作為矢量，r的方向是垂直進(jìn)入紙面的，但我們現(xiàn)在只關(guān)注他的大小。

以下使我們將會用到的Params 對象以及參數(shù)：

26_2.png

圖26.2：一卷廁紙的示意圖，標(biāo)注了力、杠桿臂、以及產(chǎn)生的扭矩。

params = Params(Rmin = 0.02 * m,
                Rmax = 0.055 * m,
                Mcore = 15e-3 * kg,
                Mroll = 215e-3 * kg,
                L = 47 * m,
                tension = 2e-4 * N,
                t_end = 180 * s)

上面出現(xiàn)的Rmin是最小半徑，Rmax是最大半徑，L是紙的長度，Mcore是紙輥中心卷心筒的質(zhì)量，Mroll是紙的質(zhì)量，tension是小貓所施加的力，單位為N，這里選擇了一個產(chǎn)生合理結(jié)果的值。

在 http://modsimpy.com/moment你可以找到簡單幾何體的轉(zhuǎn)動慣量，我將紙輥中心卷心筒視為一個“圓柱筒薄殼”，將紙輥視為一個“兩端開口的厚壁圓柱管”。

薄殼的轉(zhuǎn)動慣量是mr2 ,，其中m是質(zhì)量，r是殼的半徑。

厚壁管的轉(zhuǎn)動慣量為

I = \frac {\pi \rho h} 2 ({r_2}^4 - {r_1}^4)

式中ρ為物料密度，h為筒高，r2為外徑，r1為內(nèi)徑。

因為紙卷的外徑會隨著小貓的鋪開而變化，我們必須計算每個時間點的轉(zhuǎn)動慣量，它是一個以當(dāng)前半徑r為自變量的函數(shù)。下面是實現(xiàn)這個功能的函數(shù):

def moment_of_inertia(r, system):
    Mcore, Rmin = system.Mcore, system.Rmin
    rho_h = system.rho_h
    
    Icore = Mcore * Rmin**2
    Iroll = pi * rho_h / 2 * (r**4 - Rmin**4)
    return Icore + Iroll

rho_h是密度和高度的乘積，ρh是單位面積上的質(zhì)量。rho_h在make_system中計算:

def make_system(params):
    L, Rmax, Rmin = params.L, params.Rmax, params.Rmin
    Mroll = params.Mroll
    
    init = State(theta = 0 * radian,
                 omega = 0 * radian/s,
                 y = L)
    area = pi * (Rmax**2 - Rmin**2)
    rho_h = Mroll / area
    k = (Rmax**2 - Rmin**2) / 2 / L / radian
    
    return System(params, init=init, area=area,
                  rho_h=rho_h, k=k)

make_system同樣使用公式24.4計算k.

26_3.png

圖26.3：一個悠悠球的示意圖，標(biāo)注了重力與繩子產(chǎn)生的張力、張力的杠桿臂、以及產(chǎn)生的扭矩。

在本書的資料中，你會找到一個名為kitten.ipynb的筆記，其中包含本案例的初始代碼。用它來實現(xiàn)這個模型，并檢查結(jié)果是否可信。

26.6 模擬一個悠悠球

假設(shè)你拿著一個溜溜球，有一段繩子纏繞在它的軸上，你放下它，同時保持繩子的一端不動。當(dāng)重力使溜溜球向下加速時，繩上的張力產(chǎn)生向上的力。由于這個力作用于一個與質(zhì)心偏移的點上，它產(chǎn)生一個使溜溜球旋轉(zhuǎn)的力矩。

圖26.3是溜溜球上的力和由它們產(chǎn)生的扭矩的示意圖。外面的陰影區(qū)域是溜溜球的球身。里面的陰影區(qū)域是卷起來的弦，它的半徑隨著溜溜球展開而變化

在這個模型中，我們不能獨立計算出直線加速度和角加速度;我們需要解一個方程組:

\sum F = ma
\sum \tau = I \alpha

這里的和表示我們在把力和力矩相加起來。

在前面的例子中，線速度和角速度是相關(guān)的，因為弦展開的方式:

\frac {dy} {dt} = - r \frac {d\theta}{dt}

在這個例子中，線加速度和角加速度有相反的正負(fù)號，當(dāng)悠悠球逆時針旋轉(zhuǎn)時，θ會增加而y，也就是繩子滾動的部分的長度減小了。

對兩邊求導(dǎo)得到線加速度和角加速度的類似關(guān)系:

\frac {d^2 y} {d t^2} = -r \frac {d^2 \theta} {d t^2}

我們可以更簡潔地寫為：

a = -r \alpha

這種關(guān)系不是普遍適用的自然規(guī)律；它只適用于這種特定情況，即一個物體沿著另一個物體滾動而不拉伸或滑動。

由于我們設(shè)置問題的方式，y實際上有兩種含義:它表示滾動的弦的長度和溜溜球的高度，溜溜球的高度隨著溜溜球的下落而減小。同理，a表示滾弦長度和溜溜球高度的加速度。

我們可以通過將線性力相加來計算溜溜球的加速度:

\sum F = T - m g = m a

T是正的，因為拉力向上，而mg是負(fù)的，因為重力向下。

因為重力作用在質(zhì)心上，它不產(chǎn)生扭矩，所以唯一的扭矩是由張力產(chǎn)生的:

\sum \tau = T r = I \alpha

正(向上)的張力產(chǎn)生正(逆時針)的角加速度。

現(xiàn)在我們有三個方程，三個未知數(shù)，T, a，和α，和I, m, g，和r作為已知的量。用手解這些方程是很簡單的，但是我們也可以使用SymPy求解:

T, a, alpha, I, m, g, r = symbols('T a alpha I m g r')
eq1 = Eq(a, -r * alpha)
eq2 = Eq(T - m*g, m * a)
eq3 = Eq(T * r, I * alpha)
soln = solve([eq1, eq2, eq3], [T, a, alpha])

結(jié)果為：

T = m g I / I^*

a = mgr^2 / I^*

α = m g r / I^*

I?是增大的轉(zhuǎn)動慣量，I+mr2。要模擬這個系統(tǒng)，我們不需要T;我們可以把a(bǔ)和α直接代入到斜率函數(shù)中。

在本書的資料中，你會找到一個名為yoyo.ipynb的筆記，其中包含這些方程的推導(dǎo)過程好本案例的初始代碼。用它來實現(xiàn)并測試這個模型。

本書的中文翻譯由南開大學(xué)醫(yī)學(xué)院智能醫(yī)學(xué)工程專業(yè)2018級、2019級的師生完成，方便后續(xù)學(xué)生學(xué)習(xí)《Python仿真建?！氛n程。翻譯人員（排名不分前后）：薛淏源、金鈺、張雯、張瑩睿、趙子雨、李翀、慕振墺、許靖云、李文碩、尹瀛寰、沈紀(jì)辰、迪力木拉、樊旭波、商嘉文、趙旭、連煦、楊永新、樊一諾、劉志鑫、彭子豪、馬碧婷、吳曉玲、常智星、陳俊帆、高勝寒、韓志恒、劉天翔、張藝瀟、劉暢。

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

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

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

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容