## 11.1 計算數學 計算數學是關于通過計算來解決數學問題的科學。這里所說的“計算”既包括數值計算， 也包括符號計算；這里所說的“數學問題”可能來自純數學，更可能是從各個科學和工程領 域抽象出來的。計算數學包括很多分支，其中最核心、應用最廣的是數值方法。 數值方法 數值方法（numerical method，也稱計算方法、數值分析等）是利用計算機進行數值計 算來解決數學問題的方法，其研究內容包括數值方法的理論、分析、構造及算法等。很多科 學與工程問題都可歸結為數學問題，而數值方法對很多基本數學問題建立了數值計算的解決 辦法，例如線性代數方程組的求解、多項式插值、微積分和常微分方程的數值解法等等。 數值方法的構造和分析主要借助于數學推導，這是數學思維占主導的部分。例如，一元 二次方程的求根公式實際上給出了方程的數值解法，該公式完全是通過數學推導得出的；而 通過對該公式的分析，可以了解實數根是否存在等情形。如果問題不存在有限的求解代數式， 可以通過數學推導來尋求能得到近似解的代數式，例如將積分轉化為求和。 數值方法最終要在計算機上實現，這是計算思維占主導的部分。有人也許會認為，對于 數值計算問題，只要有了求解問題的數學公式，再將這些公式翻譯成計算機程序，問題就迎 刃而解，所以數值方法的關鍵是數學推導，而計算思維在其中并沒有什么作用。是不是這樣 呢？仍以一元二次方程 ax2+bx+c=0 的求解問題為例。這個問題的求解求根公式是已知的：  這個公式可以直接了當地翻譯成 Python 程序（程序 3.5）： ``` import math a, b, c = input("Enter the coefficients (a, b, c): ") discRoot = math.sqrt(b * b - 4 * a * c) root1 = (-b + discRoot) / (2 * a) root2 = (-b - discRoot) / (2 * a) print "The solutions are:", root1, root2 ``` 下面是此程序的一次執行結果： ``` Enter the coefficients (a, b, c): 1,-(9+10**18),9*10**18 The solutions are: 1e+18 0.0 ``` 可見，計算機求解方程 x\*\*2 -(9+10\*\*18)x + 9x10\*\*18 = 0 所給出的根是 10\*\*18 和 0，而非正確的 10\*\*18 和 9。對于這個結果，傳統的數學是無法解釋的，只有明白了計算機的能力和限制，才能給出解釋。計算思維在計算方法中的意義，由此可見一斑。 利用數值方法解決科學與工程問題大體要經過三個步驟。第一步是為問題建立數學模型，即用合適的數學工具（如方程、函數、微積分式等）來表示問題；第二步是為所建立的 數學模型選擇合適的數值計算方法；第三步是設計算法并編程實現，這里要著重考慮計算精 度和計算量等因素，以使計算機能夠高效、準確地求解問題。在計算機上執行程序得到計算 結果后，若結果不理想，多半是因為所選數值方法不合適，當然也可能是數學模型不合適。 在模型正確、編程正確的前提下，計算結果完全取決于數值方法的選擇。 本節只簡單介紹計算機的能力和限制是如何影響計算方法的選擇的。 誤差 正如前述一元二次方程求解例子所顯示的，一個正確的數學公式在計算機上卻得不到正 確的、精確的結果，這種現象主要是由誤差引起的。科學與工程計算中的誤差有多種來源， 其中建立數學模型和原始數據觀測兩方面的誤差與計算方法沒有關系，與計算方法有關的是 截斷誤差和舍入誤差。 截斷誤差是在以有限代替無限的過程中產生的，例如計算 ex 的泰勒展開式  時只能選取前面有限的 n 項，得到的是 ex 的近似值，前 n 項之后的部分就是截斷誤差。 舍入誤差是因計算機內部數的表示的限制而導致的誤差。在計算機中能夠表示的數與數 學中的數其實是不一樣的：計算機只能表示有限的、離散的數，而數學中的數是無限的、連 續的。以有限表示無限，以離散表示連續，難免造成誤差。例如 Python 中有如下出人意料 的數值計算結果： ``` >>> 1.2 - 1 0.19999999999999996 ``` 由于浮點數內部表示的限制，1.2 - 1 的結果并非精確的 0.2。又如，積分計算問題  是連續系統問題，由于計算機不能直接處理連續量，因此需要將連續的問題轉化為離散的問 題來求解。一般常用離散的求和過程來近似求解積分①。  舍入誤差的控制 計算機內部對數的表示構成一個離散的、有限的數集，而且這個數集對加減乘除四則運算是不封閉的，即兩個數進行運算后結果會超出計算機數集的范圍。這時只好用最接近的數 來表示，這就帶來了舍入誤差。因此，應當控制四則運算的過程，盡量減小誤差的影響。 在加減法運算中，存在所謂“大數吃小數”的現象，即數量級相差較大的兩個數相加減 時，較小數的有效數字會失去，導致結果中好像沒做加減一樣。例如 ``` >>> 10**18 + 9.0 1e+18 ``` > ① 據說積分號就是從 S（sum）演變而來的符號。 由此可知，當有多個浮點數相加減時，應當盡量使大小相近的數進行運算，以避免大數 “吃”小數。例如，設 x1 = 0.5055×104，x2 = x3 = ... = x11 = 0.4500（假設計算機只能支持 4 位有效數字），要計算諸 xi 的總和。一種算法是將 x1 逐步與 x2 等相加，這樣每次加法都是大 數加小數，按計算機浮點計算的規則：x1 + x2 ＝ 0.5055×104 + 0.000045×104 ＝ 0.505545×104＝0.5055×104，即產生了舍入誤差 0.45。如此執行 10 次加法之后，結果仍然是 0.5055×104，誤差積累至 10×0.45 ＝ 4.5。另一種算法是讓相近數進行運算，如 x11 + x10 = 0.9000， 在一直加到 x1，執行 10 次加法之后得到總和 0.5060×104，沒有舍入誤差。這個例子再次顯 示了“次序”在計算中的重要意義：數學上毫無差別的兩種次序在計算機中卻帶來截然不同 的結果，就像我們在第 3 章中計算 231-1 時采用 230-1+230 這個次序一樣。 當兩個相近的數相減時，會引起有效數字的位數大大減少，誤差增大。為了避免這種結 果，通常可以改變計算方法，將算式轉化成等價的另一個計算公式。例如：  在除法運算中，應當避免除數接近于零，或者除數的絕對值遠遠小于被除數的絕對值的 情形，因為這兩種情形都會使舍入誤差增大，甚至使結果溢出。解決辦法仍然是轉化為等價 算式。例如：  這里，不同計算公式的選擇就如同上述不同計算次序的選擇，雖然在數學上結果是一樣 的，但在計算機中卻存在很大差別。 計算量 站在計算機的角度，對數值方法主要關注的是算法的效率和精度。算法的效率由算法復 雜度決定，數值方法中通常用浮點乘除運算（flop）的次數來度量算法效率，稱為算法的計 算量。計算量越小，效率就越高。 當一個算法的計算量很大，并不意味著它能提高計算結果的準確度，相反倒有可能使舍 入誤差積累得更多，可謂費力不討好。利用數學推導來簡化計算公式，或者利用計算機的運 算及存貯能力來巧妙安排計算步驟，都可以減少計算量，使計算更快、更準確。 例如，設 A、B、C 分別是 10×20、20×50、50×1 的矩陣，我們來考慮如何計算 ABC。 一種算法是先算 AB，再乘 C，計算量為 10500flops；另一種算法是先算 BC，再用 A 乘， 計算量為 1200flops。顯然后一種算法大大優于前一算法，再次顯示了“次序”的妙處。 又如，考慮如何計算 x64。一種算法是將 64 個 x 逐步相乘，計算量為 63flops；另一算 法利用  其中 x\*\*2k（k＝2，4，8，16）的計算都可以利用前一步算出的結果，即  這樣計算量可以降至 10flops。 有些數值算法甚至會使計算量大到失去實際意義的地步，就如 Hanoi 塔問題的算法對較 大問題規模不可行一樣。例如求解 n 元線性方程組的克萊默法則對于較大 n 就是不可行的方 法，因為其計算量是(n+1)(n-1)(n!)+n；而高斯消去法的計算量僅為 n3/3+n2-n/3，是非常高 效的算法。 病態與良態問題 有些問題的解對初始數據非常敏感，數據的微小變化會導致計算結果的劇烈變化，這種 問題稱為病態問題，反之稱為良態問題。例如多項式 p(x) = x2+x-1150 在 100/3 和 33 處的值 分別為-5.6 和-28，數據變化只有 1%，而結果變化了 400%。又如下面這個方程組  的解是 x1 ＝ x2 ＝ x3 ＝1，當將各個系數舍入成兩位有效數字，與原來的系數雖然差別不 大，但方程組的解卻變成了 x1 ≈ -6.22，x2 ＝38.25，x3 ＝ -33.65。 相反，下面這個方程組  的解為 x1 =2，x2 = -2。若對其常數項-2 做微小擾動改為-2.005，則解變成 1.999 和-2.002， 與原來的解差別很小。可見這個問題是良態的。 數值方法主要研究良態問題的數值解法。由于實際問題的數據往往是近似值，或者是經 過舍入處理的，這相當于對原始數據的擾動，如果求解的是病態問題，則會導致很隱蔽的錯 誤結果。病態問題在函數計算、方程求根、方程組求解中都存在，它的計算或求解應當使用 專門的方法，或者轉化為良態問題來解決。 數值穩定性 求解一個問題的數值方法往往涉及大量運算，每一步運算一般都會產生舍入誤差，前面 運算的誤差也可能影響后面的運算。一個數值方法如果在計算過程中能將舍入誤差控制在一 定范圍內，就稱為數值穩定的，否則稱為數值不穩定的。例如，考慮下面這個積分的計算：  根據上面這個遞推式，可得出迭代算法：   這個算法是不穩定的，因為 I0 的舍入誤差會隨著迭代過程不斷傳播、放大。編程計算一下 可見，結果中甚至出現了負數，而根據原積分式可知 In 應該總是大于 0。 ``` >>> def f(): x = 0.1823 print "I0 =",x for n in range(1,101): x = -5 * x + 1.0 / n print "I"+str(n)+" =",x >>> f() I0 = 0.1823 I1 = 0.0885 I2 = 0.0575 I3 = 0.0458333333333 ... I97 = 1.36042495942e+63 I98 = -6.80212479709e+63 I99 = 3.40106239854e+64 I100 = -1.70053119927e+65 ``` 現在利用下列關系式  先對足夠大的 n 取 I_n 的估計值，然后再計算 I_n-1、I_n-2、…、I₁。迭代算法如下：  這個算法可使誤差逐漸減小，因此是數值穩定的，下面程序的運行結果驗證了這一點。此例又一次顯示了次序的重要性。 ``` >>> def g(): x = 0.001815 print "I100 =",x for n in range(100,0,-1): x = -x/5 + 1.0/(5*n) print "I"+str(n-1)+" =",x >>> g() I100 = 0.001815 I99 = 0.001637 I98 = 0.0016928020202 I97 = 0.00170225592249 ... I3 = 0.043138734089 I2 = 0.0580389198489 I1 = 0.0883922160302 I0 = 0.182321556794 ``` 綜上所述，數值方法以利用計算機進行數值計算的方式來解決科學和工程中抽象出來的 數學問題。與純數學方法不同，數值計算方法的構造和算法實現必須考慮計算機的能力和限 制，亦即計算思維的原則對計算方法具有重要影響。