Python科學計算（原書第2版）
Python for Scientists, Second Edition
[英] 約翰·M. 斯圖爾特（John M. Stewart） 著
江紅 餘青松 譯
|第1章|
導 論
本書英文書名為“Python for Scientists”，其具體含義是什麼呢？字典将“Python”定義為：①一種來自亞洲或者撒哈拉非洲的無毒蛇；②一種計算機腳本語言。很顯然，本書針對的是第二種定義（第二種定義的具體含義将在稍後詳細闡述）。這裡的“科學家”（Scientists），是指任何使用定量模型通過處理預先收集的實驗資料來獲得結論，或者基于一個更抽象的理論來對可能觀察到的結果進行模組化的人員，并且他們會提出“假設分析”。假設使用不同的方式來分析資料，結果如何？假設改變了模型，結果如何？是以，此處的“科學家”還包括經濟學家、工程師、數學家，以及其他通常概念下的科學家。考慮到潛在的資料量或者諸多理論模型的複雜性（非線性），使用計算機來回答上述問題正迅速成為必需。
計算機硬體的進步意味着可以以越來越快的速度來處理海量的資料和越來越複雜的模型。這些進步也意味着成本的降低，進而使得今天幾乎每個科學家都能夠使用個人計算機，即桌面工作站或者筆記本電腦，并且這兩者之間的差別正在迅速縮小。似乎也可以假定存在合适的軟體，使得這些“假設分析”問題能夠很容易地得到解答。然而，事實并非總是如此。一個明顯且實際的原因在于，雖然硬體改進有着巨大的市場，但是科學家隻占很小的一部分，是以幾乎沒有資金激勵來改進科學軟體。但對于科學家來說，這個問題很重要，我們需要更詳細地研究它。
1.1 科學計算軟體
在我們具體讨論可用的科學計算軟體之前，首先需要強調的是所有計算機軟體都分兩種類型：商用軟體和開源軟體。商用軟體是由商業公司提供的。商業公司需要同時支付工資和稅收，并為他們的股東提供投資回報。是以，商業軟體産品需要收取費用，并且為了保護其資産免受競争對手的侵權，他們不會告知客戶其軟體的工作原理。是以，最終使用者幾乎沒有機會去改變或者優化産品以滿足自己的使用要求。由于工資和稅收是經常性支出，是以公司需要經常釋出軟體更新和改進的收費通知（丹麥稅賦效應）。開源軟體則免費或者隻收取少量費用（用于媒體、郵資等）。開源軟體通常是由計算機行家個人開發的，他們經常為大學或者類似的組織工作，為同僚提供服務。開源軟體基于反版權許可進行釋出。反版權許可不給予任何人版權，也不允許将其用于商業利益。傳統經濟學可能認為開源軟體的範疇應該低于其對應的商用軟體，否則商業公司就會失去市場。然而我們将發現，事實上并非如此。
接下來，我們需要區分兩種不同類型的科學計算軟體。計算機是按照極其有限并且晦澀難懂的指令進行操作的。程式設計語言在某種程度上是人類語言的有限子集，通常是由人類編寫的指令序列，由計算機進行閱讀和了解。最常用的語言能夠表達非常複雜的數學概念，但是學習曲線非常陡峭。隻有少數語言家族（如C語言和Fortran語言）被廣泛接受，但是它們有許多不同的實作版本，如Fortran77、Fortran90、Ansi C、C++等。編譯器将人類編寫的代碼翻譯成機器代碼，進而對速度進行優化，然後進行處理。是以，它們相當于F1賽車。F1賽車中的“佼佼者”具有驚人的速度表現，但駕駛它們并不容易，需要大量的訓練和經驗。注意，編譯器需要軟體庫來補充，這些軟體庫實作了常用的數值處理算法，以及使用者通常所需要的圖形包。快速通用的庫程式包通常是昂貴的，盡管良好的公共領域程式包已經開始出現。
通常情況下，賽車并不是去超市的最好選擇，因為速度并不是最重要的因素。同樣，編譯語言并不總是适合嘗試新的數學思想。是以，對于計劃學習本書的讀者而言，除非有強制要求，否則直接使用編譯器可能沒有吸引力。是以，我們來考察另一種類型的軟體，通常稱為“科學計算軟體包”。商用軟體包包括Mathematica和Matlab，與其相對應的開源軟體包括Maxima、Octave、R和SciLab。它們都以類似的方式運作。每種軟體都提供一種特有的程式設計語言，可在使用者界面上完成輸入。在輸入了一組連貫的語句（通常隻是一個單獨的語句）之後，軟體包自動生成等效的核心語言代碼并進行動态編譯。是以，可以立即向使用者報告錯誤或結果。這種類型的軟體包被稱為“解釋器”，老讀者也許（可能帶有複雜的感覺）會想起BASIC語言。對于小型項目，與全編譯的代碼相比其慢速操作方式會被目前微處理器的速度所掩蓋，然而在處理大型任務時，其速度劣勢則會變得更明顯。
這些軟體包具有吸引力的原因至少有兩點。首先是處理資料的能力。例如，假設x是一個實變量，并且存在一個（可能未知的）函數y(x)。同時假設對于x的離散執行個體的有序集X，我們已經計算出相應的y的執行個體的集Y。随後，使用類似于plot(X, Y)的指令，該指令會立即在螢幕上顯示出完美的格式化圖形。事實上，由Matlab生成的圖形甚至可以達到出版物的品質要求。第二個優點是具有一些商用軟體顯而易見的功能，包括代數和分析功能，并且還可以将這些功能與數值和圖形功能結合起來。這些軟體包的缺點是指令語言的古怪文法和有限的表達能力。與編譯語言不同，使用這些軟體包進行程式設計通常十分困難，因為軟體包的作者并沒有考慮這些。
最好的商業軟體包非常易于使用，具有豐富的線上幫助和富有條理的文檔，這些都是對應的開源軟體所無法比拟的。然而，商業軟體包的一個主要缺點是其許可證所收取的高昂價格。大多數商業軟體包會以較低的價格提供“學生版”以便鼓勵學生熟悉該軟體包（但隻有在那些學生接受全日制教育期間才可使用）。其實這種慷慨是由其他使用者承擔的。
讓我們總結一下自身所處的環境。一方面，我們有傳統的用于數值處理的編譯語言，這些語言非常通用、非常快速，但非常難學，并且不易與圖形或者代數過程互動。另一方面，我們有标準的科學計算軟體包，這些軟體包雖然擅長将數值處理、代數和圖形互相內建，但速度較慢，并且範圍有限。
一個理想的科學計算軟體包應該具備哪些特性呢？可能包括：
1.一種既易于了解又具有廣泛表達能力的成熟程式設計語言；
2.有機內建代數、數值處理與圖形函數等功能；
3.能夠生成以速度優先運作的數值算法，其速度為編譯語言生成的數值算法中速度最快的數值算法的數量級；
4.具有一個使用者界面以提供豐富的線上幫助和完善的文檔說明；
5.相關教科書中具有豐富的擴充内容，以幫助有求知欲的讀者深入了解其概念；
6.開源軟體，免費提供；
7.在所有的标準平台（如Linux、Mac OS X、UNIX、Windows）上實作；
8.簡潔的軟體包，即便在較低配置的硬體上也可以實作。
不幸的是，迄今為止，還沒有一個“科學計算軟體包”能很好地滿足以上所有标準。
例如，如果考慮代數處理功能要求，那麼目前有兩個成熟的開源軟體包值得考慮，分别是wx-Maxima和Reduce，二者都具有強大的代數處理功能。然而，Reduce不滿足以上的第4個标準，并且二者都不滿足第3個和第5個标準。不過，在經驗豐富的使用者手中，它們都是非常強大的工具。通過附加的SymPy庫（請參見第7章），Python幾乎達到了代數能力的高标準。SageMath滿足了上面列出的除最後一條外的所有标準，它完全基于Python及其附加元件，同時還包括wx-Maxima。有關詳細資訊，請參閱第7章。是以，合理的政策是先掌握Python。如果Python中存在的少量弱點對讀者的工作來說是至關重要的，那麼就去研究SageMath。絕大多數科學家将在Python語言中找到很多有用的功能。
1991年，吉多·範羅蘇姆（Guido van Rossum）建立了Python，作為一個開源的、與平台無關的通用程式設計語言。從本質上而言，Python是一種非常簡單的程式設計語言，包括大量的附加子產品庫，以及對底層作業系統的完全通路。這意味着它可以管理和操作基于其他（甚至編譯的）軟體包建構的程式，即它是一種腳本語言。這種多功能性導緻其被超級使用者（如谷歌）以及真正的開發者隊伍采用。這也意味着它對科學家而言是一個非常強大的工具。當然，還有其他腳本語言，例如JavaTM和Perl?，但是它們的多功能性和使用者基礎都沒有滿足上述第3～5個标準。
十年前，還不可能推薦使用Python進行科學計算工作。開發人員隊伍的龐大規模意味着存在用于數值和科學應用的幾個互不相容的附加元件。幸運的是，理性占了上風，現在隻有一個數值處理附加元件NumPy，以及一個科學計算附加元件SciPy，開發人員圍繞這兩個元件完成了整合。當編寫本書的第1版時，Python中用于代數處理的SymPy庫正處于快速開發階段，是以并沒有将SymPy庫的相關知識包含在書中。雖然SymPy還沒有完全達到wx-Maxima和Reduce的能力，但它現在已經可以可靠地處理許多代數任務。
1.2 本書的規劃
本書特意寫得短小精悍，旨在向科學家們展示使用Python實作和測試相對複雜的數學算法的便捷性。我們特意側重于采用簡潔的方法來全面覆寫相關知識，目的是讓好奇的讀者盡快上手。我們的目标是為讀者建立一個處理許多基本（其實并不簡單的）任務的完善架構。顯然，大多數讀者都需要針對其特定的研究需求來進一步深入研究這些技術。但閱讀完本書後，讀者應該能打下一個堅實的基礎。
本章和附錄A将讨論如何建構Python科學計算環境。原始的Python解釋器非常基礎，其替代品IPython則更加易于使用，并且具有更加強大和多樣化的功能。第2章将緻力于采用動手操作的實踐方法來讨論IPython。
本書後續章節的内容如下。在介紹每個新特性時，我們首先試圖通過簡單的基本示例來描述，并在适當情況下通過更深入的問題來闡述。作者并不能提前預見每一個潛在讀者的數學知識水準，但在後續章節中，我們将假定讀者對基本微積分（例如一階泰勒級數）有所了解。然而，對于這些擴充問題，我們将簡述了解它們所需的背景知識，并為進一步閱讀提供适當的參考。
第3章對核心Python語言中科學家可能最感興趣的那些方面進行了簡單但相當全面的概述。Python是一種面向對象的語言，自然适合面向對象的程式設計（OOP），但這可能是大多數科學家所不熟悉的。我們将稍微涉及面向對象程式設計的主題。我們将指出3.5節中引入的容器對象在C語言或者Fortran語言中并沒有精确的類比。同樣，在3.9節中對Python類的簡要介紹可能對這兩個語言族的使用者來說并不熟悉。此章最後給出了埃拉托色尼篩選法（Sieve of Eratosthenes算法）的兩個實作，這是一個經典的問題：枚舉所有小于給定整數n的素數。一個簡單的實作需要17行代碼，但是一旦n＞105，執行時間就會非常長。然而，稍作思考并且使用已經描述的Python特性，便可實作一個更短的13行代碼的程式，它的運作速度快3000倍。但是一旦n＞108，就會耗盡記憶體（在我的筆記本電腦上）。這個練習的重點在于，選擇正确的方法（Python經常提供很多方法）是Python數值計算成功的關鍵。
第4章将通過附加子產品NumPy來擴充核心Python語言，以實作對實數和複數的高效處理。在底層使用隐藏的C/C++過程，并以接近編譯語言的速度來執行重複任務。其重點是使用向量化的代碼，而不是基于傳統的for循環或者do循環結構。向量化代碼聽起來很複雜，但正如我們将要展示的，它比傳統的基于循環的方法更容易編寫。這裡我們也将讨論資料的輸入和輸出。首先，我們讨論NumPy如何讀取和寫入文本檔案（适合人工閱讀的資料）及二進制資料。其次，我們略微涉及資料分析。我們還将總結各種函數，并簡要介紹Python的線性代數功能。最後，我們更簡要地概述另一個附加子產品SciPy，它極大地擴充了NumPy的應用範圍。
第5章将介紹附加子產品Matplotlib。該子產品受Matlab軟體包中引人注目的圖形功能的啟發，目的是對二維(x, y)繪圖進行仿真或者改進。事實上，本書後續章節中的圖形基本上都是使用Matplotlib繪制的。通過一系列的示例來說明其功能之後，我們将使用一個擴充的示例來結束此章，該示例包含功能齊全的49行代碼，用于計算并生成曼德爾布羅特集（Mandelbrot set）的高清晰度繪圖。
第6章将讨論擴充到三維圖形的難點，例如幾何面z=z(x, y)的繪制。從某種程度而言，可以使用Matplotlib子產品來處理，但對于更一般的情況，則需要調用Mayavi附加子產品，此章将簡單介紹該子產品，同時給出一些示例代碼。如果使用這樣的圖形是讀者的主要興趣所在，那麼需要進一步研究這些子產品。
最後一個介紹性章節（第7章）将展示SymPy的代數功能，盡管它存在局限性，但是讀者可能會收獲驚喜。
如果讀者已經具備一些Python的經驗知識，那麼當然可以跳過這些章節的部分内容。然而，這些章節的主旨在于實際操作的方法。是以強烈鼓勵讀者嘗試運作相關的代碼片段。一旦讀者了解了這些代碼，便可以通過修改它們來加深了解。這些“黑客”實驗（“hacking” experiment）取代了傳統教科書中所包含的練習。
上述章節涵蓋了Python為增強科學家的計算機經驗而提供的基本工具。接下來的章節将包含哪些内容呢？
上述章節的明顯遺漏之處在于沒有涵蓋廣泛的資料分析主題（僅在4.5節簡短闡述）。這主要有如下三個原因：
1.最近出現了一個叫作Pandas的附加子產品，該子產品使用NumPy和Matplotlib來解決資料分析的問題，而且它帶有詳細且全面的文檔（這在4.5節中進行了描述）。
2.Pandas子產品的作者之一寫了一本書（McKinney（2012）），該書概述了IPython、NumPy和Matplotlib，并詳細闡述了Pandas應用程式。
3.這不是我的工作領域，是以隻會改述源代碼。
是以，我選擇專注于科學家的模組化活動。一種方法是隻針對生物資訊學、宇宙學、晶體學、工程學、流行病學或金融數學等某一門學科中的問題。事實上，本書的前半部分内容可以衍生出一系列書籍，如“面向生物資訊學的Python”等。另一種不那麼奢侈但更有效的方法（即本書後半部分采用的方法）是針對上述所有領域（以及更多領域）。這依賴于數學的統一性：在某個領域中的問題一旦簡化到核心無量綱形式後，往往看起來與另一個領域中簡化到無量綱形式的問題相類似。
這種特性可以通過下面的例子來說明。在種群動态學中，我們可能研究一個單物種，它的種群數量N(T)依賴于時間T。給定充足的食物供應，我們可能期望指數增長dN/dT=kN(T)，其中生長常量k具有維數1/T。然而，這種增長通常存在限制。包含這些的簡單模型是“邏輯斯谛方程”（logistic equation）。
它允許一個穩定的常數種群N(T)=N0。許多教科書中對這一方程的生物學背景進行了讨論，如Murray（2002）。
在（同質球型對稱）宇宙學中，密度參數Ω依賴于尺度因子a：
其中w通常假設為常量。
目前數學生物學和宇宙學沒有太多的共同之處，但很容易看出式（1-1）和式（1-2）代表相同的方程。假設我們使用比例t=kN0T來縮放方程（1-1）中的自變量T，進而得到一個新的無量綱的時間坐标t。同樣，我們引入無量綱變量x=N/N0，進而使得方程（1-1）成為邏輯斯谛方程。
在一般的相對性理論中，沒有任何理由去選擇或者偏向某一種時間坐标。是以，我們可以選擇一個新的時間坐标t：a=et/(1+3w)。然後設定x=Ω，我們将看到方程（1-2）也簡化為方程（1-3）。是以，相同的方程可以在許多不同的領域出現。在第8～10章中，為了簡潔明了，我們将使用類似方程（1-3）的極小方程。如果讀者所要解決的問題的最簡形式與上述代碼片段類似，則很自然可以套用這段代碼，以處理求解的問題所對應的原本冗長的代碼。
第8章将讨論涉及常微分方程的四類問題。我們首先簡要介紹解決初值問題的技術，然後給出若幹示例，包括兩個經典的非線性問題：範德波爾振蕩器和洛倫茲方程。接下來，我們将考察兩點邊值問題，并研究線性Sturm-Liouville特征值問題，以及一個有關非線性Bratu問題的拓展練習。延遲微分方程在控制理論和數學生物學中經常出現，例如邏輯斯谛方程和麥克-格拉斯（Mackey-Glass）方程，是以下一節将讨論它們的數值求解方法。在這一章的最後，我們将簡要介紹随機微積分方程和随機常微分方程。特别是，我們将讨論一個與金融數學中的布萊克-斯克爾斯（Black-Scholes）方程緊密相關的簡單例子。
這裡還将介紹兩個與科學家相關的重要Python主題。第一個主題是如何嵌入使用其他語言編寫的代碼，這包括兩個方面的内容：①重用現有的遺留代碼（通常用Fortran編寫）；②如果性能分析器表明程式代碼執行速度受若幹Python函數的影響而嚴重減慢，如何重新使用Fortran或者C語言編寫并替換這些“惹是生非”的函數？第二個主題是科學計算使用者如何能夠有效使用Python的面向對象程式設計（OOP）特性？
第9章通過一個擴充的例子來解決第一個主題。我們首先了解如何使用僞譜方法來解決由偏微分方程支配的大量演化問題，即初始值或者初始邊值問題。為了簡潔起見，我們隻讨論涉及一個時間次元和一個空間次元的問題。這裡我們将闡述周期性空間相關性的問題可以使用傅裡葉方法非常有效地進行處理，但對于更一般的問題，則需要使用切比雪夫變換。然而在這種情況下，還沒有令人滿意的Python黑盒可用。而事實上，必要的工具已經在傳統的Fortran77代碼中編寫實作。附錄B中列出了這些功能，并且我們将展示如何用極少的Fortran77知識來構造速度極快的Python函數以完成所需的任務。我們的方法依賴于NumPy f2py工具，而它包含在所有推薦的Python釋出版本中。如果讀者對重用先前存在的遺留代碼感興趣，那麼即使這一章處理的示例與讀者面臨的問題無關，也值得研究。有關f2py的其他用途，請參閱1.3節。
從科學家的角度來看，面向對象程式設計最有用的特征之一是類的概念。類存在于C++（而不是C）和Fortran90及其後版本（而不是Fortran77）中。然而，這兩種語言的實作都十分複雜，是以通常初級程式員往往會回避使用。相比之下，Python的實作要簡單得多，而且更加友好，代價是省略了其他語言實作的一些更神秘的特性。我們在3.9節中對其文法進行了簡要介紹。然而在第10章，我們将給出一個更加實用的例子：使用多重網格（multigrid）來求解任意維數的橢圓型偏微分方程，盡管為了簡明起見，示例代碼是針對二維網格的。多重網格目前是一個經典問題，其最佳描述是使用遞歸方式進行定義，是以我們将使用若幹篇幅來描述它，至少能概述清楚。先前存在的遺留代碼相當複雜，因為作者需要用不支援遞歸的語言（例如，Fortran77）來模拟遞歸。當然，我們可以使用第9章中闡述的f2py工具來實作這個代碼。此外，我們也可以使用Python類和遞歸來構造一個簡單的多重網格代碼。作為一個具體執行個體，我們使用了Press等（2007）中對應章節内的樣本問題，以便感興趣的讀者可以比較非遞歸方法和面向對象的方法。如果讀者對多重網格并沒有特别的興趣，但确實關注涉及關聯數學結構的問題，并且這些問題經常出現在諸如生物資訊學、化學、流行病學和固态實體學等領域，那麼理所當然應該仔細閱讀最後這一章，以了解如何在數學上精确地解釋問題，進而輕松地構造Python代碼來求解問題。
1.3 Python能與編譯語言競争嗎
對Python及其科學計算軟體包最常見的批評是它們在處理複雜的實際問題時與編譯代碼相比速度太慢。注重運作速度的讀者可能需要了解最近的一項研究，即通過各種方法來應對簡單的“數字處理”問題。雖然其最後的結果圖是在單個處理器上處理一個特定問題，但它們确實給出了性能上“大緻正确”的印象。他們使用完全編譯的C++程式來解決這個問題，并将其速度作為一個基準。使用第3章所述技術（即核心Python技術）的解決方案速度慢了約700倍。一旦讀者使用浮點子產品NumPy和第4章中描述的技術，代碼速度便會慢大約10倍左右，并且與Matlab的性能估計差不多一緻。然而，正如研究指出的，有很多方法可以将Python加速到C++性能的80%左右。其中一些實踐在計算機科學中卓有成效。
特别是，存在一個對科學家非常有用的工具：f2py。我們将在第9章中詳細讨論f2py，闡述如何重用傳統的Fortran遺留代碼。它還可以用于通路标準的Fortran庫，例如NAG庫。另一個用途是加速NumPy代碼，進而提高性能。為了了解其工作原理，假設我們已經開發了一個程式（例如本書後面幾節中描述的那些程式），程式使用了大量的函數，每個函數執行一個簡單的任務。這個程式運作正常，但速度慢得令人無法接受。注意，可以直接擷取Python代碼所需的詳細的時間性能資料。Python包括一個可以在工作程式上運作的“性能分析器”（profiler）。“性能分析器”按執行函數所花費的時間順序來輸出函數的詳細清單。“性能分析器”的使用非常便捷，具體請參見2.5節。通常，會存在一到兩個函數，它們耗費很長的時間來執行簡單的算法。
這就是f2py的切入點。因為函數很簡單，即使初學者也可以很快建立等價語言（例如，Fortran77或者Ansi C）的代碼。此外，因為我們所編碼的内容十分簡單，是以不需要對應語言（例如，Fortran95或者C++）的進階功能（這些需要花費精力去學習）。接下來，我們使用f2py工具将代碼封裝在Python函數中，并将它們嵌入Python程式内。隻要有一點經驗，我們就可以達到與完全使用其他語言（例如，Fortran95）編寫的程式相同的運作速度。
1.4 本書的局限性
全面闡述Python及其各個分支需要大量的篇幅，而且在成書之前也許就已經過時了。是以，本書的目标是為讀者提供一個起點，使讀者掌握如何使用基本的附加軟體包。一旦掌握了一些有關Python應用的經驗知識，就可以進一步探索感興趣的領域了。
我自己也意識到本書未能涉及一些非常重要的概念，例如，雙曲問題的有限體積方法、并行程式設計和實時圖形等，Python可以有效地應用于這些領域。在科學研究的前沿有一支非常強大的Python開發者隊伍，通過網際網路很容易通路他們的工作成果。請讀者把本書看作通往科學研究前線的交通工具吧。
1.5 安裝Python和附加軟體包
Matlab和Mathematica的使用者習慣于使用定制的內建開發環境（IDE）。在啟動螢幕中，使用者可以使用内置編輯器審閱代碼，編寫、編輯和儲存代碼片段，并運作具體的程式。由于作業系統Mac OS X和大多數Linux版本事實上都包含核心的Python版本，是以許多計算機人員和其他經驗豐富的黑客會告訴使用者，安裝額外的軟體包是一件很簡單的事情，而且使用者可以在短短一小時内開始編碼，進而彌補了一些相對于IDE的劣勢。
不幸的是，專家們的觀點是錯誤的。本書中涉及的Python系統會将語言運作到極限，是以所有的附加軟體包必須彼此相容。與許多其他作者一樣，我本人也經曆了數小時的挫折來嘗試專家們的政策。請讀者閱讀附錄A，以節省時間和精力。基于上述原因，附錄A主要針對初學者。
當然，這裡也包含一定的争議和困擾（盡管不多且不太嚴重）。那麼權衡之計是什麼呢？如果讀者遵循附錄A中建議的路線，則最終會擁有一個無縫工作的系統。由于原始的Python解釋器确實不那麼友好，是以所有的IDE提供商都提供了一個“Python模式”，但是他們聲稱所提供的需求已經被增強的解釋器IPython所取代。實際上，在其最新版本中，IPython希望超越Matlab、Mathematica，以及商業IDE中與Python相關的功能。特别是，其允許使用者使用自己喜歡的編輯器，而不是他們的編輯器，并且根據使用者的需要定制指令。附錄A和第2章中将闡述這些内容。