我们已经看到了原理上可以处理等离子体问题的一种方法，其基本思路就是用网格上的有限差分法来求解电磁场；用集团代表点来计算带电粒子运动；用蒙特卡洛抽样方法来处理所有量子力学相关过程。那么在知道这些东西之后，剩下的就是把方法实现为程序，然后拿去处理实际问题。

在这里你会立刻意识到一个问题，这个程序的开发模型是自顶向下的，所以传统软件工程的那些古老模型和管理方案都可以直接使用。事实上差不多就是如此，你几乎可以按三十年前软件工程教科书上的方法来建立工作计划，而且(!)可以按照“人月妄想”【注1】来管理进度。

当然你会问，既然这问题已经到了这种完全是管理艺术的层次，还要做什么科研？正统的回答是这样：难道因为操作系统设计都已经工业化了，操作系统的研究就应该被列入“毫无创新性”的黑名单吗？当然实际上远远没有这么简单。问题的关键之一，是上面所说的几个部分，其实全都是尚未解决--而且也可能长期无法解决的问题。

我们先说一个不那么关键但需要理解的事情：这个方法的几个部件并不是精密耦合的，你完全可以在理论允许的情况下换掉上面的一个模块。最典型的情况，带电粒子之间的相互作用原则上是电磁场，应该用麦克斯韦方程组描述，但是假如问题涉及到的频率比较低，也不存在明显的电磁波辐射出来，那你把电磁场求解换成静电场的泊松方程甚至库仑定律也可以处理问题【注2】。粒子的运动方程本来应该是相对论的洛伦兹方程式，但如果涉及到的能量很低，你直接用牛顿方程就简单很多。

当然这一切都依赖于你准备拿这个程序做什么。脱离应用目标（或者说你的市场所在）来讨论算法，在其他领域也许是有好处的，但在我们这个特定领域，几乎无法实现【注3】。所以下面给大家一个例子，看看一个比较典型的问题的开发过程。【注4】

问题本身非常简单：假设一个航天器在空间中飞行，太阳风吹拂出的电子和离子轰击在航天器的表面上，我们通常认为电子和离子击中物质表面就会附着在上面。因为通常来说电子比离子飞行速度快得多，在密度相同的情况下，单位时间落在表面上的电子电量就要大于离子电量，于是航天器主要充的是负电。随着充电的进行，航天器带负电压会排斥电子，于是电子充电量会降低。等到电子和离子充电流量相等的时候，航天器电压趋于稳定，最后形成一个稳态的负电压。这个电压会和电子平均动能一个量级，所以正常情况下，深空或者高地球轨道航天器会带一个很高的负电压。由于遮挡问题，不同位置的电压有所不同，在有些情况下，这种电压差很大，会损坏航天器的内部电子电路。我们需要一个程序完成两件事情：1：根据设备电压反推背景环境电子和离子能量；2：在已知电子和离子背景参数的时候求出航天器的最高电压差。

暂时不提这问题的一些特殊细节，显然我们首先需要一个能解算等离子体运动的基本程序，或者说求解器。我们可以估计出航天器的充电电流是非常小的，时间演变也非常缓慢（譬如说，从毫秒到天的特征时间），电子能量也就是几千到几万电子伏特，所以可以认为整个问题就是一个静电问题。下面我们简单的来考虑一下问题。

既然用粒子云网格方法，那么网格肯定是必要的，问题在于航天器本身是个实体，所以我们需要把这个实体分解到网格上。最简单的选择就是空间分成一系列网格，每个网格给个描述，是实体还是空间。好吧当然，这分割够粗糙，问题主要是实体外形是个曲面，我们简单的用一堆小方块去替换它，那肯定会有精度问题。然而这是程序细化所要解决的问题。

接下来我们需要计算带电粒子的运动，这无非就是求解每个粒子的牛顿运动方程：dV/dt=F/m。这看上去只是个简单的常微分方程数值解的问题，虽然涉及到的方程多了一点（通常我们用至少一千万个模拟粒子），但反正方程极其简单。事实上当然不会如此简单，因为这里处理的问题包含静电简谐振荡，而我在很早就说过，对于最简单的简谐振荡问题，大家用了几十年的4阶龙格-库塔方法本质上是错误的：它会带来能量不守恒问题。这问题本身是NASA的“计算专家”熟知的，并且后来被专业数学家发现然后建立了一套称为“辛算法”的技巧。这里有一个好消息和一个坏消息。坏消息是辛算法的思路在现在这个问题上过于复杂了，因为这不完全是一个常微分方程组问题：我们同时还要求解电磁场方程，而电磁场总是个偏微分方程，所以你需要完全从头推导辛算法那些理论（大概这五年才有人做成这件事）；好消息是在很多情况下，我们不是很关心格式的阶次，因此可以退回到2阶；2阶的时候最简单的“蛙跳格式”就是保证能量守恒的。所以我们首先可以用蛙跳格式来处理问题。

现在你看到这个问题一个巨大的（也是根本性的简化）简化：求解带电粒子运动之后，我们才来求解电磁场方程，严格的说这么做是不自洽的：带电粒子和电磁场应该耦合在一起求解，只不过那么做要产生一个维度高得难以理解的非线性微分方程（你问这东西是椭圆还是抛物还是双曲的？别逗了，这是无限维相空间中的方程，鬼才知道它的具体特性），没人知道要怎么求解。现在的做法，电磁场方程已经退化成了一个标准的静电场Poisson方程，于是可以套用能在网格上求解Poisson方程的任何方法：有限差分，有限元，甚至拟谱展开，快速多极子什么的。（啊，有限元！你一定会立刻想到，这么成熟的方法，看来问题已经简单了，根本不需要再絮絮叨叨了。。。。）

下面仍然是一个坏消息：你可以去找一本偏微分方程数值求解的专著，然后会发现大量篇幅在讲一个问题如何离散成线性方程组，这里不准备细说，总之你这么理解，反正是个微分方程，一个地方的空间导数和函数值的关系，这就在各个格点的函数值之间建立了关系，最后形成一个N维线性方程组，N等于网格格点数。通常发明一种离散方法框架就可以写教科书的一章，问题在于，这种线性方程组怎么解？

当然你会说，直接matlab就好。对于教科书上的典型问题，比如结构力学有限元之类，确实如此，所以他们通常都说，这种计算和写程序的问题是小trick，只有“物理和数学”才是重要的。但这里有个很特殊的地方：从理论上我们会知道，PIC算法每个时间步应该不高于等离子体振荡周期的三分之一，这可能对应也就是几个纳秒或者皮秒的时间，而你要算到微秒甚至毫秒级别，也就是你要求解微分方程几万到几百万次。而问题的规模呢？就算用最粗糙的网格，我们也需要至少每个方向100个网格，于是对应网格数是一百万个。好吧，一百万个未知量的线性代数求解要多久？（你可以用matlab求解一个一万乘一万的线性方程组，差不多1秒？因此100万的大概要一万秒）不妨假定你能在100秒内求解出来（试试看，这不是个很简单的事情），然后一百万步就是一亿秒，三年多。

所以你看到问题就这么简单，算的够不够快的问题。事实上线性方程组的求解算是基础算术中最成熟的问题了，最基本的数学书会教你克莱姆法则（行列式那个法则），当然实际上根本没人能用那个法则，那个法则的计算量是未知数个数N的阶乘数，这个概念你自己去想；其次是高斯消元法，大概是沈括时代发明的东西，计算次数大约是未知数个数的三次方，听起来已经是“多项式级别”的，不过实际上还是不太现实，毕竟一百万的三次方是1e18，还是太扯了。由于这个方程组是“稀疏的”，有一些专门的技巧来进行逼近，称为迭代算法，在一般的“数值线性代数”或者“应用线性代数”书上会给你介绍常见的迭代法，比如最传统的高斯-赛德尔迭代，或者“比较近代化”的共轭梯度(CG)算法，甚至GMRES 之类。不过实际上你会发现，在其他方面无往而不利的CG/GMRES 算法，在这里基本上是废的：首先对于这么大的网格经常不收敛；其次就算收敛也会慢到无法接受（比如可能要迭代几百次），也许那些比较“工程化”的学科可以忍受几百次迭代，但是在这里是完全不可想象的。

当然，问题并不是无解，在很久以前（大概80年代）有几个“数学家”（我不确定真正的数学家是不是把他们当数学家，毕竟还有认为数学物理都不算数学或者在数学中是不重要的外道的牛人）发明了一种方法，简单的说就是认为Poisson方程这种东西，解应该是“逐步细化的”，比如说你可以先用粗糙的网格估计一下结果，然后再计算细网格应该修正多少。这办法能出奇迹：对于Poisson问题，这种方法能在很多情况下实现最多10次迭代就必然给出结果。

当然，实践中这方法有一系列问题，最大问题是，稍微改变一下边界条件或者方程结构，这办法就不行了。但是架不住“伪数学家”能凑合，我们可以把这个方法当成一个预估手段应用到共轭梯度方法上，这就是MG-Precond-CG/GMRES 方法。这种方法兼顾了可靠性和速度，算是一种强有力的方法。

遗憾的是，很多时候你得自己写这种代码。偏偏写这类代码还不是一个很小的工作，没有几万行代码是下不来的。所以这就成了国家实验室和高校的锅。SANDIA之类首当其冲，他们为了干自己的活，专门纠集一帮人开发这种算法，当然这算法开发完了就对他们没用了，所以他们把这些东西放在网上，当成一种开放程序库给研究者提供方便（很遗憾国内的大学和研究所还没有这种习惯）。

当然，这种方法是一种比较通用的方法，通用就意味着不够专门，比如说，同样一百万的未知量，用这些人写的通用求解器需要上百秒钟，而专门针对特定的背景，可以先把方程化简，那也许只要零点一秒。所以实际中你必须自己权衡，究竟是先写完基本框架，然后用人家写好的方程求解器补进去，还是自己一直开发下去。权衡的核心是，直接用通用求解法能不能满足需求，以及自己开发一个需要多少工作量（包括是否要克服某些数学困难）【注5】。在现在的情况下，一百秒钟的求解速度是很难满足要求的，不过我们还有变通。

总之，现在基本问题已经清楚了，你要做出框架，然后对几个部分调研找出可行的解决方法。这里的可行性是主要的问题，因为很多文献中或者教科书中“可行”的方法其实是不具备可行性的。如果你依靠经验确认所有的可行性都不是问题，那么就可以开始下一步的工作；如果你发现有些东西可行性不足（比如说，常规的线性方程组求解法对于现在这个问题就存在可行性问题），那你需要开发一组专门的算法来处理。

回到这个问题，我们看到这个问题主要是巨大的计算量，导致事情不那么简单。从而在可行性论证这个阶段，我们就必须去把问题的难点扫除。由于数学上已经没有什么技巧可用，我们只能再下沉到计算机技术这个层面。

【注1】《人月神话》，一本很有意思的书。

【注2】库伦定律和麦克斯韦方程组之间的关系很有意思。

【注3】这种专门应用和通用理论之间的关系，常常是科学家和工程师以及教科书作者之间互相鄙视的主要话题。

【注4】这是我们组的近期工作之一，基本上是在NSFC课题的基础上发展出来的应用。我们的NSFC课题本身看上去和这个问题毫无关系（一个关于放电过程的模拟课题），但这个问题涉及到的各种技术几乎都是来源于NSFC课题的论文。

【注5】估计工作量的唯一方法是做一个相关的问题确认一下，而通用求解法能不能解决问题更麻烦，它需要(1)理论上排除可能的“不适用性”(2)实践中做一个小规模的同类问题测试一下。所以对于高校教师和研究生，这都是一种考验。