计算机科学丛书
点击查看第二章
点击查看第三章
Python科学计算（原书第2版）
Python for Scientists, Second Edition
[英] 约翰·M. 斯图尔特（John M. Stewart） 著
江红 余青松 译
|第1章|
Python for Scientists, Second Edition
导 论
本书英文书名为“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章中将阐述这些内容。