中學(xué)時(shí)，沒(méi)有什么比數(shù)學(xué)圖紙更恐怖的了。

許多問(wèn)題沒(méi)有現(xiàn)成解析方案，還有一些問(wèn)題雖然可以解決，但需要大量的計(jì)算工作。這種情況下，你需要算法來(lái)輔助。

希望你沒(méi)被算法整吐。萬(wàn)一你真的吐了，也要看到好的一面：又可以再吃頓午餐了:]

在這個(gè)教程中，你將學(xué)習(xí)到算法的概念，以及如何使用它們來(lái)解決難以分析的問(wèn)題。通過(guò)playground，很容易使解決方案可視化，易于查看。

即使你不是個(gè)數(shù)學(xué)愛(ài)好者，對(duì)物理或計(jì)算機(jī)科學(xué)也不大感興趣，你仍會(huì)在這個(gè)教程中找到一些有意思的東西。你只需要有一些微積分和基礎(chǔ)物理學(xué)的知識(shí)。

你將學(xué)到如何使用數(shù)值算法解決兩個(gè)困難問(wèn)題。但是為了講的更清楚點(diǎn)，這兩個(gè)問(wèn)題也可通過(guò)解析法解決。算法在解析法無(wú)法工作時(shí)更為理想，即便如此，通過(guò)比較兩種方法，也能更加容易的理解它們的本質(zhì)。

數(shù)值算法是什么？

簡(jiǎn)單說(shuō)來(lái)，數(shù)值算法是一類(lèi)解決數(shù)學(xué)問(wèn)題的方法，它們不依賴于閉式解析解。

問(wèn)題來(lái)了，什么是閉式解析解？

若有一個(gè)公式，我們可以使用已知數(shù)，通過(guò)有限的操作步驟，可以獲得一個(gè)準(zhǔn)確的結(jié)果，即稱之為閉式解析解。

再簡(jiǎn)單點(diǎn)來(lái)理解一下：如果你可以使用代數(shù)的方式，找到一個(gè)表達(dá)式來(lái)解決一個(gè)未知量問(wèn)題，代入某些已知數(shù)即可得到結(jié)果，就說(shuō)明你已經(jīng)得到了一個(gè)閉式解析方法。

何時(shí)使用數(shù)值算法？

許多問(wèn)題沒(méi)有現(xiàn)成解析方案，還有一些問(wèn)題雖然可以解決，但需要大量的計(jì)算工作。這種情況下，你需要算法來(lái)輔助。

例如，假定你需要編寫(xiě)一個(gè)物理引擎，用來(lái)計(jì)算許多物體在有限時(shí)間內(nèi)的行為。這時(shí)你就可以使用數(shù)值算法來(lái)更快地得到結(jié)果。

這樣做也有副作用，更快的計(jì)算速度意味著結(jié)果不會(huì)十分精確。然而，大多數(shù)情況下，這樣的結(jié)果已經(jīng)夠用了。

天氣預(yù)報(bào)就得益于數(shù)值計(jì)算。天氣變化的速度、影響因素的數(shù)量都是十分驚人的。這是一個(gè)非常復(fù)雜的系統(tǒng)，只有數(shù)值模擬才能完成預(yù)知未來(lái)的任務(wù)。

可能正是因?yàn)槿狈@些算法，你的iPhone總是告訴你要下雨了，而外面還是艷陽(yáng)高照。

[[142780]]

開(kāi)始

作為熱身活動(dòng)，我們來(lái)玩?zhèn)€游戲，接下來(lái)，你將計(jì)算出一個(gè)給定數(shù)字的平方根。這兩個(gè)方法都將使用二分法（bisection method）。神奇的是，你可能已經(jīng)知道了這個(gè)方法，但是不知道它的名字。

回想一下，在兒童時(shí)代，我們可能玩過(guò)這樣的游戲：選中100以內(nèi)的一個(gè)數(shù)字，另外一個(gè)人來(lái)猜。你只會(huì)提示他猜的數(shù)字大了或者小了。

游戲開(kāi)始。小明告訴小強(qiáng)開(kāi)始猜，小強(qiáng)說(shuō)我猜1，小明說(shuō)小了，小強(qiáng)又猜2，小明說(shuō)小了。小強(qiáng)接下來(lái)再猜3,4...最后選中了5，終于對(duì)了。

5步就猜中了，不錯(cuò)。但是如果小明選的是78，猜中就需要花點(diǎn)時(shí)間。

這個(gè)游戲如果使用二分法（bisection method）來(lái)玩的話，猜中的速度會(huì)快很多。

二分法

我們知道數(shù)字介于1和100之間，因此除了一個(gè)一個(gè)的猜，或者隨便瞎猜外。我們把數(shù)字分為兩個(gè)區(qū)間：a = [1,50], b = [51, 100]。

接下來(lái)我們判斷數(shù)字是介于哪一個(gè)區(qū)間，這很簡(jiǎn)單，只用問(wèn)數(shù)字是不是50。如果數(shù)字小于50，那么區(qū)間b就不用考慮了。接下來(lái)我們繼續(xù)把a(bǔ)區(qū)間分成兩半，再重復(fù)這個(gè)步驟。

例如：

數(shù)字是60.區(qū)間為：a = [1,50], b = [51, 100]。

第一步，小強(qiáng)說(shuō)50（也即是a區(qū)間的上限），小明說(shuō)小了。這時(shí)小強(qiáng)就知道了數(shù)字肯定在b區(qū)間上。

第二步，分解b區(qū)間為： c=[51,75]，d=[76,100]。還是選擇c區(qū)間的上限75，小強(qiáng)告訴他大了。這就意味著數(shù)字肯定在c區(qū)間上。

繼續(xù)分解。。。

使用這個(gè)方法，7次嘗試即可得到正確答案，一個(gè)一個(gè)試則需要60次。

1. 50 -> 小了

2. 75 -> 大了

3. 62 -> 大了

4. 56 -> 小了

5. 59 -> 小了

6. 61 -> 大了

7. 60 -> 對(duì)了！

計(jì)算x的平方根，過(guò)程也類(lèi)似。數(shù)字x的平方根介于0和x之間。也可以記做：(0,x]。如果數(shù)字x大于或等于1，可以記做：[1, x]。

分解區(qū)間，得到a=(0, x/2]，b=(x/2, x]。

如果數(shù)字x是9，區(qū)間是[1, 9]，分解后的區(qū)間為a=[1, 5]，b=(5, 9]，中間值m為(1+9)/2 = 5。

下一步，檢查m * m - x，是否大于期望的精度。在這里，我們檢查m * m 大于或小于等于x。如果大于，我們使用(0, x/2]區(qū)間繼續(xù)檢查，否則，使用(x/2, x]區(qū)間。

看一下實(shí)際步驟，初始化m=5, 期望準(zhǔn)確值為0.1:

1. 計(jì)算m * m - x: 5 * 5 - 9 = 11

2. 檢查結(jié)果是否小于等于期望值：25 - 11 <= 0.1?

3. 顯然不滿足，繼續(xù)下一個(gè)區(qū)間：m * m是否大于9？是。接下來(lái)使用區(qū)間[1, 5]，測(cè)試值為(1+5)/2=3。

4. 計(jì)算m * m - x：３* 3 - 9 = 0

5. 查查是否滿足期望值：9 - 9 <= 0.1?

6. 搞定。

備注：你可能想知道括號(hào)是什么意思？區(qū)間有固定的格式（下限和上限）。不同的括號(hào)代表不同的區(qū)間邊界。圓括號(hào)表示邊界不在區(qū)間范圍內(nèi)（即開(kāi)區(qū)間），而方括號(hào)表示邊界在區(qū)間范圍內(nèi)（閉區(qū)間）。如(0, a] 包含a而不包含0.在上面的例子中：區(qū)間 (0, a/2]包含a/2而不包含0；區(qū)間(a/2, a]表示所有大于a/2, 小于a的數(shù)。

在Playground上使用二分法

現(xiàn)在，是時(shí)候應(yīng)用學(xué)到的理論了。我們來(lái)自己實(shí)現(xiàn)二分算法。創(chuàng)建一個(gè)新的playground文件，添加如下的代碼：

func bisection(x: Double) -> Double { 
//1 
    var lower = 1.0 
    //2 
    var upper = x 
    //3 
    var m = (lower + upper) / 2 
    var epsilon = 1e-10 
    //4 
    while (fabs(m * m - x) > epsilon) { 
    //5 
        m = (lower + upper) / 2 
        if m * m > x { 
            upper = m 
        } else { 
            lower = m 
        } 
    } 
   
    return m 
} 

1. 設(shè)置區(qū)間下限為1

2. 設(shè)置區(qū)間上限為x

3. 找到中間值，并定義期望精確度

4. 檢查操作是否能滿足精確度

5. 如果不滿足，找到新的中間值，定義新的區(qū)間，繼續(xù)查找。

添加如下代碼來(lái)測(cè)試該函數(shù)：

let bis = bisection(2.5) 

在 m = (lower + upper) / 2這一行的右邊，可以看到代碼執(zhí)行了35次，意味著找到結(jié)果需要35步。

#p#

馬上我們就能看到playground一個(gè)可愛(ài)功能帶來(lái)的好處：數(shù)值的完整歷史都可以查看。

二分法可以成功的算出實(shí)際結(jié)果的近似值。通過(guò)數(shù)值的歷史記錄圖，就可以看到數(shù)值算法是如何逼近正確的解。

按下option+cmd+enter打開(kāi)輔助編輯器，點(diǎn)擊m = (lower + upper) / 2代碼行右邊的圓按鈕在輔助編輯器中添加歷史記錄。

你會(huì)看到方法一點(diǎn)點(diǎn)的跳轉(zhuǎn)到正確結(jié)果上。

古典數(shù)學(xué)也搞的定

下一個(gè)算法需要追溯到古代，它起源于公元前1750年的古巴比倫，在公元100年前亞歷山大的希羅所著《度量論》（Metrica ）一書(shū)中有所描述。這就是為何它被稱作“希羅方法”。可以通過(guò)這個(gè) 鏈接 了解更多關(guān)于希羅的知識(shí)。

這個(gè)方法使用函數(shù)，這里你想要計(jì)算a的平方根。如果你能找到函數(shù)曲線的“0值點(diǎn)（zero point）”，這個(gè)點(diǎn)上的函數(shù)值為0，那么求出x值就能得到a的平方根。

要完成這項(xiàng)工作，我們首先選擇任意x值作為起始值，計(jì)算該值對(duì)應(yīng)點(diǎn)的tangent值（函數(shù)的切線），然后查看tangent線上的0點(diǎn)（即該直線和x軸的交點(diǎn)）。然后我們使用這個(gè)0點(diǎn)再次作為起始值，重復(fù)多次直到滿足精度。

因?yàn)槊坑?jì)算一次tangent值，都會(huì)更加逼近真正的0值，

這個(gè)過(guò)程會(huì)逐漸逼近真正的結(jié)果。下圖展示了求解a=9時(shí)，a的平方根，且起始值為1.

起始點(diǎn)x0 = 1，生成一條紅色的tangent線，接著生成了x1點(diǎn)，又生成了紫色的線；又生成了x2點(diǎn)，連接了藍(lán)色的線，最終找到了正確答案。

當(dāng)古典數(shù)學(xué)邂逅了Playground

我們有很多古巴比倫人沒(méi)有的好東西，比如：playground。在playground上添加如下代碼，并查看：

func heron(x: Double) -> Double { 
  //1 
  var xOld = 0.0 
  var xNew = (x + 1.0) / 2.0 
  var epsilon = 1e-10 
   
  //2 
  while (fabs(xNew - xOld) > epsilon) { 
    //3 
    xOld = xNew 
    xNew = (xOld + x / xOld) / 2 
  } 
   
  return xNew 
} 

這段代碼做了什么？

1. 創(chuàng)建一個(gè)變量來(lái)存儲(chǔ)當(dāng)前結(jié)果。xOld是上次計(jì)算的結(jié)果，而xNew是實(shí)際結(jié)果。

• 找到初始化xNew的方法是一個(gè)良好的起點(diǎn)

• epsilon是期望的精確度

• 這個(gè)例子中，我們計(jì)算平方根精確度為小數(shù)點(diǎn)后10位。

2. 使用while循環(huán)查找是否達(dá)到期望的精確度。

3. 如果達(dá)不到精確度，設(shè)置xNew為xOld，繼續(xù)下一次迭代。

使用如下代碼驗(yàn)證該函數(shù)的作用：

     
let her = heron(2.5) 

希羅方法只需要5次迭代即可找到正確結(jié)果。

在代碼行xNew = (xOld + x/xOld)/2右邊點(diǎn)擊圓角button，添加值歷史，就能看到第一次迭代就能找到一個(gè)不錯(cuò)的近似值。

模擬諧振子運(yùn)動(dòng)

在這個(gè)章節(jié)中，我們來(lái)看看如何使用數(shù)值積分算法來(lái)模擬簡(jiǎn)諧系統(tǒng)運(yùn)動(dòng)——一種基本的動(dòng)態(tài)系統(tǒng)。

這個(gè)系統(tǒng)可以描述很多現(xiàn)象，比如鐘擺的擺動(dòng)、彈簧的振動(dòng)。特別說(shuō)來(lái)，它可以描述某段時(shí)間內(nèi)物體移動(dòng)了一定偏移量的場(chǎng)景。

這個(gè)例子中，假定有一個(gè)有質(zhì)量的物體連接在彈簧上。為了簡(jiǎn)化問(wèn)題，我們忽略掉阻力和重力。因此唯一的作用力就是彈簧的彈力，它將物體拉回到原來(lái)的位置。

在這樣的假定下，你只需要使用兩個(gè)力學(xué)定律來(lái)處理問(wèn)題：

• 牛頓第二運(yùn)動(dòng)定律, 它描述了物體上的作用力和加速度的關(guān)系。

• 胡克定律，它描述了彈力和物體偏移量之間的比例關(guān)系。這里k是彈力系數(shù)而x是物體的偏移量。

簡(jiǎn)振方程

因?yàn)閺椓κ俏矬w上唯一的作用力，我們將上面的兩個(gè)方程式整理：

簡(jiǎn)化后寫(xiě)作：

也被記作，也即是共振頻率的平方。

更精確的方程式如下：

其中A是振幅，這里即是物體的偏移量，被稱為相位差。這兩個(gè)值初始化時(shí)都被設(shè)置為定值。

如果設(shè)置時(shí)間t=0，則，且，你可以簡(jiǎn)單的計(jì)算出振幅和相位差：

來(lái)看一個(gè)例子，我們有一個(gè)質(zhì)量為2kg的物體連接在彈簧上，彈簧的彈力系數(shù)為k=196N/m。在初始時(shí)間t=0時(shí)，彈簧移動(dòng)了0.1米。我們可以通過(guò)公式計(jì)算振幅、相差和共振頻率。

在Playground上實(shí)驗(yàn)一下：

使用該公式計(jì)算任意時(shí)間點(diǎn)對(duì)應(yīng)的值之前，我們需要添加一些代碼。

回到playground文件，在最后添加如下代碼：

//1 
typealias Solver = (Double, Double, Double, Double, Double) -> Void 
   
//2 
struct HarmonicOscillator { 
  var kSpring = 0.0 
  var mass = 0.0 
  var phase = 0.0 
  var amplitude = 0.0 
  var deltaT = 0.0 
   
  init(kSpring: Double, mass: Double, phase: Double, amplitude: Double, deltaT: Double) { 
    self.kSpring = kSpring 
    self.mass = mass 
    self.phase = phase 
    self.amplitude = amplitude 
    self.deltaT = deltaT 
  } 
   
  //3 
  func solveUsingSolver(solver: Solver) { 
    solver(kSpring, mass, phase, amplitude, deltaT) 
  } 
} 

這些代碼塊中做了什么？

1. 定義一個(gè)函數(shù)類(lèi)型的的typealias，函數(shù)有5個(gè)Double類(lèi)型的參數(shù)，返回為空。

2. 創(chuàng)建一個(gè)struct來(lái)定義諧振子。

3. 這個(gè)函數(shù)只是簡(jiǎn)單的創(chuàng)建用來(lái)解振動(dòng)問(wèn)題的Clousure。（并無(wú)真實(shí)的計(jì)算代碼）

再來(lái)看一下精確方案

代碼如下：

func solveExact(amplitude: Double, phase: Double, kSpring: Double, mass: Double, t: Double) { 
  var x = 0.0 
  //1 
  let omega = sqrt(kSpring / mass) 
   
  var i = 0.0 
   
  while i < 100.0 { 
//2 
    x =  amplitude * sin(omega * i + phase) 
    i += t 
  } 
} 

這個(gè)方法包含了所有解決運(yùn)動(dòng)方程需要的參數(shù)。

1. 計(jì)算共振頻率

2. 在while循環(huán)中計(jì)算物體當(dāng)前的位置，i用作下次計(jì)算的自增變量。

添加如下測(cè)試代碼：

let osci = HarmonicOscillator(kSpring: 0.5, mass: 10, phase: 10, amplitude: 50, deltaT: 0.1) 
osci.solveUsingSolver(solveExact) 

這個(gè)方案中的方法有點(diǎn)奇怪：它有參數(shù)傳入，但不返回?cái)?shù)據(jù)，也沒(méi)有顯示任何東西。

#p#

這樣做有什么好處？

使用這個(gè)函數(shù)的目的在于，在while循環(huán)中，可以算出振動(dòng)過(guò)程中具體的動(dòng)態(tài)值。在Playground中，可以觀察這些動(dòng)態(tài)值的歷史記錄。

在x = amplitude * sin(omega * i + phase) 處添加值歷史記錄，我們就能看到運(yùn)動(dòng)的軌跡。

[[142803]]

既然第一個(gè)精確算法已經(jīng)驗(yàn)證通過(guò)，下面我們看一下數(shù)值計(jì)算方案。

歐拉方法

歐拉方法是用來(lái)求數(shù)值積分的一種方法。該方法1768年再Leonhard Euler所著Institutiones Calculi Integralis《積分學(xué)原理》一書(shū)中首次提出。要查看該方法的詳情，請(qǐng)參考：http://en.wikipedia.org/wiki/Euler_method

歐拉方法的核心思想是通過(guò)使用折線逼近曲線。

具體方法是計(jì)算一個(gè)給定點(diǎn)的斜率，然后繪制一條同樣斜率的折線。在這條折線的末尾，繼續(xù)計(jì)算斜率，繪制另外一條線。正如你所見(jiàn)，該算法的精確度取決于折線的長(zhǎng)度。

你想知道deltaT做了什么嗎？

該數(shù)值算法中需要使用一個(gè)步長(zhǎng)，這對(duì)算法的精確度很重要。大步長(zhǎng)導(dǎo)致低精確度，但是執(zhí)行速度塊。反之，精確度會(huì)提高，速度會(huì)降低。

deltaT變量就代表了步長(zhǎng)(step size)。我們初始化這個(gè)值為0.1， 意味著我們計(jì)算每0.1秒物體移動(dòng)的位置。在歐拉算法中，意味著折線在x軸上的投影長(zhǎng)度為0.1。

在編寫(xiě)代碼前，你需要再看一下這個(gè)公式：

二階微分方程可以轉(zhuǎn)化為兩個(gè)一階微分方程。可以被寫(xiě)為和

我們可以用差商來(lái)完成轉(zhuǎn)換：

也會(huì)得到：

有了這些等式，我們可以直接實(shí)現(xiàn)歐拉方法。

在solveExact方法后添加代碼：

func solveEuler(amplitude: Double, phase: Double, kSpring: Double, mass: Double, t: Double) { 
  //1 
  var x = amplitude * sin(phase) 
  let omega = sqrt(kSpring / mass) 
  var i = 0.0 
  //2 
  var v = amplitude * omega * cos(phase); 
  var vold = v 
  var xoldEuler = x 
   
  while i < 100 { 
    //3 
    v -= omega * omega  * x * t 
    //4 
    x += vold * t 
    xoldEuler = x 
   
    vold = v 
    i+=t 
    } 
} 

代碼的內(nèi)容：

1. 設(shè)置當(dāng)前的位置，和omega的值。

2. 計(jì)算當(dāng)前速度。

3. 在while循環(huán)中，通過(guò)一階函數(shù)計(jì)算出新的速度。

4. 使用速度計(jì)算新的位置，在結(jié)束處，使用步長(zhǎng)t增加i的值。

現(xiàn)在，添加以下代碼測(cè)試：

osci.solveUsingSolver(solveEuler) 

在xoldEuler = x位置添加值歷史，并查看。你會(huì)看到這個(gè)方法顯示一個(gè)振幅增加的正弦函數(shù)。因?yàn)闅W拉方法并不是一個(gè)精確算法，而這里過(guò)大的步長(zhǎng)0.1也導(dǎo)致了計(jì)算結(jié)果不精確。

[[142811]]

以下是另外一個(gè)函數(shù)圖像，步長(zhǎng)為0.01，這樣明顯更好。因此，使用歐拉方法時(shí)你要想到，步長(zhǎng)越小，結(jié)果越好。但是使用折中的步長(zhǎng)，執(zhí)行起來(lái)更為簡(jiǎn)單。

[[142812]]

速度Verlet算法(Velocity Verlet)

最后討論的算法叫做速度Verlet。它和歐拉方法的思路一樣，但是計(jì)算新位置的方式有些許差異。

歐拉在計(jì)算新位置時(shí)，忽略實(shí)際的加速度，公式為：,而速度Verlet算法在計(jì)算時(shí)考慮到了加速度：。這再步長(zhǎng)相等的時(shí)候，結(jié)果更優(yōu)。

在solveEuler方法后添加新的代碼：

func solveVerlet(amplitude: Double, phase: Double, kSpring: Double, mass: Double, t: Double) { 
  //1 
  var x = amplitude * sin(phase) 
  var xoldVerlet = x 
  let omega = sqrt(kSpring / mass) 
   
  var v = amplitude * omega * cos(phase) 
  var vold = v 
  var i = 0.0 
   
  while i < 100 { 
    //2 
    x = xoldVerlet + v * t + 0.5 * omega * omega * t * t 
    v -= omega * omega * x * t 
    xoldVerlet = x 
    i+=t 
  } 
} 

代碼的內(nèi)容：

1. 循環(huán)前的所有代碼和歐拉方法中的一樣。

2. 根據(jù)速度計(jì)算出新的位置，增加i的值。

在Playground中測(cè)試代碼:

osci.solveUsingSolver(solveVerlet) 

[[142815]]

接下來(lái)做什么？

你可以通過(guò) 此鏈接 下載完整工程.

希望你在這場(chǎng)數(shù)值世界旅行中獲得樂(lè)趣。了解一些算法，即便只是一些古典數(shù)學(xué)的趣味歷史，都會(huì)對(duì)你帶來(lái)幫助。