第二講：動力學方程的數(shù)值積分

計算機動畫的本質(zhì)上是解微分方程均践，因此本文主要介紹微分方程常用的一些解法桂躏。

particle system

首先從最簡單的內(nèi)容開始抢埋，也就是particle system宋渔。

上面是particle system的一個簡單的案例州疾，也就是一個彈簧+一個block，block在彈簧的作用下左右運動傻谁，其中彈簧不受力的情況下block所在的位置為原點孝治，這里用q表示位置，那么上述狀態(tài)可以用q=0來表示审磁，定義向左為負谈飒，向右為正，相應(yīng)的速度與加速度也可以基于下面的公式得到：

$\dot{q} = \frac{dq}{dt}$
$\ddot{q} = \frac{d\dot{q}}{dt}$

動能kinetic energy $E_k$ 與勢能potential energy $E_p$ 則由如下的公式進行表示：

$E_k = \frac{1}{2}m\dot{q}^2$
$E_p = \frac{1}{2}kq^2$

其中m與k分別表示block的質(zhì)量與彈簧的彈性系數(shù)态蒂，而彈簧的彈力可通過下述公式計算得到：
$f = -kq$

根據(jù)前面的公式可以知道杭措，彈力是conservative force的一種。

所謂的conservative force指的是將一個particle從一個點移動到另一個點做的功是跟移動的路徑無關(guān)的力钾恢，在這個定義下手素，如果將一個particle沿著任意路徑移動最后回到原點，那么力的整體做功（微分下可以用 $\int f(x) * dx$ 表示）為0瘩蚪，在conservative force的作用下泉懦，不論沿著哪條路徑，只要保證首尾點是相同的疹瘦，那么都會導(dǎo)致勢能的變化幅度是相同的崩哩，常見的力中，重力跟彈力是conservative force，而摩擦力則不是conservative force邓嘹，從這個角度來理解酣栈，conservative force是跟位置相關(guān)的力，而這種力作用的效果就是將位置的變化轉(zhuǎn)換為勢能的變化汹押。

一個力場（force field）如果滿足下面三個公式中的任意一個矿筝，那么就可以說明這個力場是conservative的：

力的curl是0向量
$\vec{\nabla} \times \vec{F} = 0$
在2D空間中，上述公式可以簡化為：
$\frac{\partial F_y}{\partial x} - \frac{\partial F_x}{\partial y} = 0$
閉環(huán)移動時棚贾，做功為0
$W = \oint \vec{F} \mathrmygz0onm \vec{r} = 0$
力可以通過對勢能的梯度取反來表示：
$\vec{F} = -\vec{\nabla}\Phi$

前面兩個公式計算比較復(fù)雜窖维，通常我們會使用第三個公式進行conservative force的判斷，根據(jù)前面彈簧粒子的公式鸟悴，我們可以判斷彈力就是conservative force陈辱。

再說回彈簧粒子系統(tǒng)，根據(jù)能量守恒定律细诸，我們有：
$\frac{1}{2}m\dot{q}^2 + \frac{1}{2}kq^2 = C$
其中C是常量，對這個公式的左右兩邊對時間進行求導(dǎo)陋守，就可以得到：
$m\dot{q}\ddot{q} + kq\dot{q}=\dot{q}(m\ddot{q}+kq) = 0$
上述這個公式有兩個解：
$\dot{q} = 0 \\ m\ddot{q} = - kq$
其中第一個解的含義是速度 $\dot{q}$ 不變震贵，也就是物件靜止，這種情況我們是不關(guān)心的水评，我們主要關(guān)心第二個解猩系，這個解我們換個寫法：
$\ddot{q} = \frac{- kq}{m}$
用中文的話翻譯一下，就是物體運動的加速度大小與物體的受力成正比中燥，跟物體的質(zhì)量成反比寇甸，這就是牛頓第二定律（作為回顧：牛頓第一定律也稱為慣性定律，可以表述為任何物體都要保持勻速直線運動或靜止狀態(tài)疗涉，直到外力迫使它改變運動狀態(tài)為止拿霉；牛頓第三定律可以表述為相互作用的兩個物體之間的作用力和反作用力總是大小相等，方向相反咱扣，作用在同一條直線上）绽淘。

下面再來看一個復(fù)雜一點的案例，block從一個變成兩個：

動能公式就變成了：
$E_k = \frac{1}{2}m_1\dot{q_1}^2 + \frac{1}{2}m_2\dot{q_2}^2$

用矩陣來表示就可以寫成：
$E_k = \frac{1}{2} \left[ \begin{matrix} \dot{q_1} \\ \dot{q_2} \end{matrix} \right]^T \left[ \begin{matrix} m_1 & 0 \\ 0 & m_2 \end{matrix} \right] \left[ \begin{matrix} \dot{q_1} \\ \dot{q_2} \end{matrix} \right]$

這里需要注意的是闹伪，我們前面說過沪铭，單個block的情況下，我們是以彈簧原始長度時block所在的位置為原點的偏瓤，而實際上我們這里選擇的坐標系是比較隨性的杀怠，而不是有固定的一套標準，而使用不同的坐標系厅克，得到的公式是不一樣的赔退，雖然他們輸出的結(jié)果是一樣的（能量不隨坐標系的變化），而這給了我們一個啟發(fā)已骇，在實際使用中要謹慎選擇坐標系离钝，不同的坐標系計算公式不一樣票编，導(dǎo)致計算效率會受影響，而這對于結(jié)果的輸出并沒有任何收益卵渴。

數(shù)值積分

重寫一下前面的牛頓第二定律的公式：

$\ddot{q} = \frac{- kq}{m}$

這是一個以q為變量的二階微分方程慧域，在物理動畫模擬中，我們最終需要的是計算各個時刻的位置浪读，也就是解出每個time step下的q昔榴，而求解微分方程的方式就是數(shù)值積分。

為了解前面這個二階微分方程碘橘，我們這里嘗試將之改寫成一個一階的常微分方程（目的是降低求解的復(fù)雜度互订？），令
$y = \left[ \begin{matrix} q \\ \dot{q} \end{matrix} \right]$
那么有
$\frac{\textyrz02lky}{\text5ckcy00t} = \left[ \begin{matrix} \dot{q} \\ \ddot{q} \end{matrix} \right] = \left[ \begin{matrix} \dot{q} \\ \frac{-kq}{m} = \frac{F(q, \dot{q})}{m} \end{matrix} \right] = \mathscr{F} (y)$
上面為了進行一般化處理痘拆，用 $F(q, \dot{q})$ 替代了 $-kq$ 仰禽，那么怎么來理解這個y以及y的微分呢？

如下圖所示纺蛆，我們以q作為橫坐標吐葵， $\dot{q}$ 作為縱坐標繪制一個坐標系（這個坐標系這里稱之為phase space，當然桥氏，這里只是一個近似温峭，真實的phase space的縱坐標是動量 $m\dot{q}$ ，不過不影響）字支，當我們?nèi)↑c $y = (0, \dot{q})$ 來看凤藏，其微分結(jié)果為 $(\dot{q}, 0)$ ，也就是一條水平方向（向右）的向量堕伪，如下面所示：

同樣揖庄，如果我們?nèi)↑c $y = (q, 0)$ 來看，其微分結(jié)果為 $(0, \frac{-kq}{m})$ 刃跛，也就是一條垂直方向（向下）的向量抠艾，如下面所示：

實際上，對單block的彈簧系統(tǒng)而言桨昙，對于每個y坐標检号，在phase space中得到的y的微分，稱為y的速度向量蛙酪，都可以得到這個向量是沿著順時針方向的齐苛，如果我們做特殊處理，取k = m桂塞，那么將所有點連起來凹蜂，我們就能得到一個正圓，如下圖所示：

不同半徑的正圓對應(yīng)于彈簧和諧運動的振幅，通過這個正圓也可以看到玛痊，在q為0的時候汰瘫，速度的絕對值最大，反之速度為0的時候擂煞，則是位置的絕對值最大的時候混弥；另外，對于任意一個y对省，其速度 $\frac{\mathrmq09sp7yy}{\mathrmlxfyymst}$ 都是這個正圓上的切線蝗拿，也就是說下一刻的位置同樣是在正圓上，如果沒有摩擦力的存在蒿涎，那么就意味著這個block將在這個圓上永無止盡的運動下去哀托。

在進行數(shù)值積分方法介紹之前，我們先來看下不同的積分方法的評價標準劳秋，總的來說仓手，對于一個數(shù)值積分方法好不好，我們通常有如下三個評判標準：

stability玻淑，穩(wěn)定性俗或，整個系統(tǒng)會隨著時間的變遷變得穩(wěn)定，不論是勻速運動岁忘，還是隨著能量的衰減逐漸趨于靜止，或者在某個有限的能量的范圍內(nèi)來回振蕩
accuracy区匠，精確度干像，指的是我們通過數(shù)值積分模擬出來的結(jié)果跟真實的精確解之間的誤差
convergence，收斂度驰弄，指的是我們在進行數(shù)值積分模擬計算的時候麻汰，模擬的time step在不斷收縮的情況下，其模擬的精度是否不斷逼近真實解戚篙。

將上面的公式做一下近似五鲫，這里用explicit euler（顯式歐拉）公式進行模擬，即 $\scr{F}(y)$ 中y取k：
$\frac{y(k+1) - y(k)}{h} \approx \frac{\textttqurm5y}{\texti2pt2adt} = \mathscr{F} (y(k))$
從而得到
$y(k+1) = h \cdot \mathscr{F}(y(k)) + y(k)$

按照上面這個公式岔擂，由于 $\scr{F}(y)$ 是在y點的切線位喂，也就是說計算得到的下一個坐標點將落在當前坐標點所規(guī)范的圓形外面，那么經(jīng)過這個公式計算出來的運動軌跡將不能保證是圓形的乱灵，而是半徑不斷擴大的螺旋曲線塑崖，也就是說，這個方法是不穩(wěn)定的痛倚，即系統(tǒng)的能量越來越大（半徑對應(yīng)著振幅）：

在顯式歐拉方法下规婆，我們看下整體方案的誤差，根據(jù)泰勒展式：
$y(t_0+h) = y(t_0) + \dot y(t_0) \cdot h + \ddot y(t_0) \cdot \frac{h^2}{2} + ...$
可以看到，顯式歐拉只使用了 $\ddot y$ 之前的部分抒蚜，也就是說是屬于一階微分方程掘鄙，其誤差自然就可以算出來等于 $O(h^2)$ 了。

如果前面不用顯式歐拉嗡髓，而是使用implicit euler（隱式歐拉）操漠，即 $\scr{F}(y)$ 中y取k+1：
$\frac{y(k+1) - y(k)}{h} \approx \frac{\textttq5jxly}{\texttmm7cqlt} = \mathscr {F}(y(k+1))$

在這種情況下，右邊是一個方程器贩，在彈簧系統(tǒng)中颅夺，這個方程是可以顯式表達的，但是在其他的系統(tǒng)尤其是復(fù)雜系統(tǒng)中蛹稍，這個方程是沒有辦法表達出來的吧黄，也就是說這個公式通過解析方法是沒有辦法計算得到的。

$y(k+1) - h\cdot \mathscr{F}(y(k+1)) = y(k)$

根據(jù)上面的公式唆姐，仿造前面顯式歐拉的邏輯拗慨，我們可以看出來，y(k)是比y(k+1)半徑要大的奉芦，也就是說赵抢，在這個系統(tǒng)中，半徑在不斷減小声功，能量在不斷衰退烦却，最終衰減到0，那么就處于靜止先巴，也就是說這種方式得到的是一個穩(wěn)定的解（當然其爵，這里需要說清楚的是，不是所有的系統(tǒng)使用隱式歐拉都是穩(wěn)定的伸蚯，只有那些能量守恒的系統(tǒng)才是穩(wěn)定的）摩渺，不過呢由于 $\scr{F}$ 沒有解析解，因此求解會很麻煩剂邮。

下面要介紹的是一個叫做symplectic euler（辛歐拉唁桩，或者半顯式歐拉火脉，半隱式歐拉方案）的數(shù)值積分方法局待，根據(jù)前面的公式亡问，我們有：
$y = \left[ \begin{matrix} q \\ \dot{q} \end{matrix} \right]$
$\dot{y} = \left[ \begin{matrix} \dot{q} \\ \ddot{q} \end{matrix} \right] = \scr{F}(q, \dot{q})$
將之轉(zhuǎn)換為微分方程：
$\frac{1}{h} \left[ \begin{matrix} q(k+1) - q(k) \\ \dot{q}(k+1) - \dot{q}(k) \end{matrix} \right] = \left[ \begin{matrix} \dot{q} \\ \frac{F}{m}(q, \dot{q}) \end{matrix} \right]$
又回到前面的問題，右邊的向量中的q是應(yīng)該取k還是k+1時刻的值瑞眼，前面我們說過龙宏，取k的值是顯式歐拉的方案，好處是數(shù)值是已知的伤疙，這個方程可以解出來的银酗，缺點則是系統(tǒng)不穩(wěn)定辆影，而取k+1的時候則是隱式歐拉的方案，好處是方案是穩(wěn)定的黍特，而缺點則是方程求解困難蛙讥，基于這一點，我們這里考慮使用一個混合的方案灭衷，即上述方程中第二個分量使用顯式歐拉方案次慢，第一個分量則使用隱式歐拉方案，如下所示：
$\frac{1}{h} \left[ \begin{matrix} q(k+1) - q(k) \\ \dot{q}(k+1) - \dot{q}(k) \end{matrix} \right] = \left[ \begin{matrix} \dot{q}(k+1) \\ \frac{F}{m}(q(k), \dot{q}(k)) \end{matrix} \right]$
根據(jù)上面的公式翔曲，我們可以先根據(jù)第二個分量的公式迫像，計算出 $\dot{q}(k+1)$ ，之后這個數(shù)值可以代入到第一個公式瞳遍，求得 $q(k+1)$ 闻妓，這里借用了顯式歐拉的好處，只是需要確定的是掠械，是否同時帶入了顯式歐拉的缺點——不穩(wěn)定由缆。

在彈簧系統(tǒng)中，上述公式可以表示為：
$\frac{1}{h} \left[ \begin{matrix} q(k+1) - q(k) \\ \dot{q}(k+1) - \dot{q}(k) \end{matrix} \right] = \left[ \begin{matrix} \dot{q}(k+1) \\ \frac{-k\cdot q(k)}{m} \end{matrix} \right]$
從而求得：
$\dot{q}(k+1) = \dot{q}(k) - \frac{h \cdot k \cdot q(k)}{m} \\ q(k+1) = q(k) + h \cdot \dot{q}(k+1) = q(k) + \dot{q}(k) \cdot h- \frac{ k \cdot q(k)}{m} h^2$
根據(jù)這公式通過不斷迭代可以得到：
$q(k+1) = q(k) - \frac{kh^2}{m}\sum_{i=0}^k q(i) +\dot{q}(0) \cdot h$
好像也不能說明啥猾蒂，寫成矩陣形式：
$\left[ \begin{matrix} q(k+1) \\ \dot{q}(k+1) \end{matrix} \right] = \left[ \begin{matrix} 1-\frac{k}{m} h^2 && h \\ -\frac{kh}{m} && 1 \end{matrix} \right] \left[ \begin{matrix} q(k) \\ \dot{q}(k) \end{matrix} \right]$
這個等式右邊的2x2的矩陣均唉，我們可以很容易計算出其行列式=1，而根據(jù)線性代數(shù)的理論肚菠，我們可以知道舔箭，矩陣的行列式恰好等于變換前后兩個形狀的面積之比，也就是說蚊逢，某個區(qū)域的數(shù)據(jù)集合經(jīng)過這個矩陣變化后得到的新的區(qū)域限嫌，兩個區(qū)域的面積之比就恰好等于矩陣的行列式，而這里行列式為1就說明經(jīng)過這個矩陣的變換时捌，數(shù)據(jù)是基本維持穩(wěn)定的，不會膨脹炉抒，也不會坍縮奢讨，而是能量守恒的（自然也是穩(wěn)定的系統(tǒng)）。

而雖然上述結(jié)論是在彈簧系統(tǒng)中推導(dǎo)出來的焰薄，但實際上這個結(jié)論的應(yīng)用面是可以推廣到一般的（當然拿诸，不是所有的系統(tǒng)都可以使用），對于系統(tǒng)內(nèi)僅包含動能和勢能的系統(tǒng)而言塞茅，在經(jīng)典力學范疇中亩码，這個結(jié)論是成立的，即通過辛歐拉算法得到的解都是穩(wěn)定的（甚至能量守恒的）野瘦。

最后編輯于：2022.04.06 10:08:44

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者

人面猴
序言：七十年代末描沟，一起剝皮案震驚了整個濱河市飒泻，隨后出現(xiàn)的幾起案子，更是在濱河造成了極大的恐慌吏廉，老刑警劉巖泞遗，帶你破解...
沈念sama閱讀 212,816評論 6贊 492
死咒
序言：濱河連續(xù)發(fā)生了三起死亡事件，死亡現(xiàn)場離奇詭異席覆，居然都是意外死亡史辙，警方通過查閱死者的電腦和手機，發(fā)現(xiàn)死者居然都...
沈念sama閱讀 90,729評論 3贊 385
救了他兩次的神仙讓他今天三更去死
文/潘曉璐我一進店門佩伤，熙熙樓的掌柜王于貴愁眉苦臉地迎上來聊倔，“玉大人，你說我怎么就攤上這事生巡“颐铮” “怎么了？”我有些...
開封第一講書人閱讀 158,300評論 0贊 348
道士緝兇錄：失蹤的賣姜人
文/不壞的土叔我叫張陵障斋，是天一觀的道長纵潦。經(jīng)常有香客問我，道長垃环，這世上最難降的妖魔是什么邀层？我笑而不...
開封第一講書人閱讀 56,780評論 1贊 285
?港島之戀（遺憾婚禮）
正文為了忘掉前任，我火速辦了婚禮遂庄，結(jié)果婚禮上寥院，老公的妹妹穿的比我還像新娘。我一直安慰自己涛目，他們只是感情好秸谢，可當我...
茶點故事閱讀 65,890評論 6贊 385
惡毒庶女頂嫁案：這布局不是一般人想出來的
文/花漫我一把揭開白布。她就那樣靜靜地躺著霹肝，像睡著了一般估蹄。火紅的嫁衣襯著肌膚如雪。梳的紋絲不亂的頭發(fā)上沫换，一...
開封第一講書人閱讀 50,084評論 1贊 291
城市分裂傳說
那天臭蚁，我揣著相機與錄音，去河邊找鬼讯赏。笑死垮兑，一個胖子當著我的面吹牛，可吹牛的內(nèi)容都是我干的漱挎。我是一名探鬼主播系枪，決...
沈念sama閱讀 39,151評論 3贊 410
雙鴛鴦連環(huán)套：你想象不到人心有多黑
文/蒼蘭香墨我猛地睜開眼，長吁一口氣：“原來是場噩夢啊……” “哼磕谅！你這毒婦竟也來了私爷？” 一聲冷哼從身側(cè)響起雾棺，我...
開封第一講書人閱讀 37,912評論 0贊 268
萬榮殺人案實錄
序言：老撾萬榮一對情侶失蹤，失蹤者是張志新（化名）和其女友劉穎当犯，沒想到半個月后垢村，有當?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體，經(jīng)...
沈念sama閱讀 44,355評論 1贊 303
?護林員之死
正文獨居荒郊野嶺守林人離奇死亡嚎卫，尸身上長有42處帶血的膿包…… 初始之章·張勛以下內(nèi)容為張勛視角年9月15日...
茶點故事閱讀 36,666評論 2贊 327
?白月光啟示錄
正文我和宋清朗相戀三年嘉栓，在試婚紗的時候發(fā)現(xiàn)自己被綠了。大學時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片拓诸。...
茶點故事閱讀 38,809評論 1贊 341
活死人
序言：一個原本活蹦亂跳的男人離奇死亡侵佃，死狀恐怖，靈堂內(nèi)的尸體忽然破棺而出奠支，到底是詐尸還是另有隱情馋辈，我是刑警寧澤，帶...
沈念sama閱讀 34,504評論 4贊 334
?日本核電站爆炸內(nèi)幕
正文年R本政府宣布倍谜，位于F島的核電站迈螟，受9級特大地震影響，放射性物質(zhì)發(fā)生泄漏尔崔。R本人自食惡果不足惜答毫，卻給世界環(huán)境...
茶點故事閱讀 40,150評論 3贊 317
男人毒藥：我在死后第九天來索命
文/蒙蒙一、第九天我趴在偏房一處隱蔽的房頂上張望季春。院中可真熱鬧洗搂，春花似錦、人聲如沸载弄。這莊子的主人今日做“春日...
開封第一講書人閱讀 30,882評論 0贊 21
一樁弒父案，背后竟有這般陰謀
文/蒼蘭香墨我抬頭看了看天上的太陽宇攻。三九已至惫叛，卻和暖如春，著一層夾襖步出監(jiān)牢的瞬間逞刷，已是汗流浹背挣棕。一陣腳步聲響...
開封第一講書人閱讀 32,121評論 1贊 267
情欲美人皮
我被黑心中介騙來泰國打工，沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留亲桥，地道東北人。一個月前我還...
沈念sama閱讀 46,628評論 2贊 362
代替公主和親
正文我出身青樓固耘，卻偏偏與公主長得像题篷，于是被迫代替她去往敵國和親。傳聞我的和親對象是個殘疾皇子厅目，可洞房花燭夜當晚...
茶點故事閱讀 43,724評論 2贊 351

第二講：動力學方程的數(shù)值積分

particle system

數(shù)值積分

推薦閱讀更多精彩內(nèi)容