1. 引言
无力场是天体物理中经常用到的磁场模型，尤其在太阳物理的有关研究中受到极大的关注 [1] 。数学上，无力场是指满足如下方程的向量场：
(1)
(2)
其中，是空间位置的标量函数，称作无力因子。当恒为零时，无力场退化为势场，即无电流场；当恒为非零常数时，称作线性无力场；当恒为非常数时，称作非线性无力场，此时无力场方程为二阶非线性偏微分方程组，只有少数情况下可以得到一些特解。1990年，Low和Lou [2] 首次给出了一类多极非线性无力场，在天体物理领域获得广泛的关注和应用，详见Wiegelmann和Sakurai的综述文章 [3] 。
在轴对称条件下，任何磁场都可以在球坐标系中写成如下形式：
(3)
其中，与是和的函数，与方位角无关。结合这一条件，Low和Lou将无力场方程转化为如下形式的常微分方程：
(4)
其中，为自变量；是的函数，也是有待求解的未知函数，满足边界条件。
方程(4)是带有未知参数的二阶非线性常微分方程，Low和Lou将其视作特征值问题进行了“数值研究”，并给出了情形下的两组解，但他们在文章中省略了具体的求解过程，由此带来一些困惑，而方程(4)本身的研究也停滞了十多年。直到最近，随着有关研究的深入，非线性无力场的Low-Lou解法才再次引起人们的兴趣 [4] - [7] 。这些进展都是在Low-Lou解法的框架内展开的，由于天体物理领域的作者们侧重于物理内容，求解过程多有省略。迄今为止，对于方程(4)的数值求解，在方法上仍然没有明朗化。本文提出一种参数打靶法，作为对Low-Lou解法的探讨和补遗，并给出更多可选的Low-Lou解。
2. 参数打靶法
考虑两点边值问题：
(5)
其中，为参数。若已知，则可以采用通常的打靶法求解该问题。首先将原方程转化为如下的初值问题：
(6)
其中，为在初始点处的斜率。则问题转化为寻找合适的，使得该初值问题的解满足原边值问题的右端边界条件。其次，设相应于的初值问题(6)的解为，并且假定关于连续，则要找的为满足如下方程：
(7)
上述方程可用牛顿迭代法求解。
从另一个角度看，数值求解方程(7)的过程，相当于在一定范围内调整，使得误差函数的值接近于零。依此观点，当为未知参数时，可以固定，然后在一定范围内调整，使得误差函数的值接近于零，由此获得原边值问题的数值解。这里的关键在于根据原始边值问题的上下文，选取合适的固定值。我们将这种方法称为参数打靶法，见算法2.1。显然，参数打靶法可看作打靶法的变体形式或灵活运用。
笔者注意到，最近杨旦旦等人 [8] 也表述了一种处理未知参数的打靶法，并用于火箭和卫星中的液体燃料晃动问题的求解。本文提出的参数打靶法，主要是针对特征值(或参数值)不唯一的情形。特别地，我们首先在打靶区间内获得完整的误差函数曲线，然后一次性求出该曲线的所有零点，可有效地防止“漏算”的情况发生。
3. 算例及分析
为了便于对照Low和Lou的原始结果，我们考虑方程(4)中的情形，即
(8)
其中，为自变量；是的函数，满足边界条件。特别地，Low和Lou将此方程作为特征值问题求解，其中看作特征函数，看作特征值。由于方程没有常数项，边界条件又是齐性的，从而特征函数允许所谓的自由规范化。固定起见，Low和Lou规定：在处，取。为了实施参数打靶法，我们需要选定参数的打靶区间，这里令。进一步，在算法2.1中，用表示，取，。通过试算可知，特征值的分布越来越稀疏。为了获得精确结果，可将打靶区间分为若干段，逐段实施参数打靶法；计算结果参见图1及表1。注意，特征值的实际分布是越来越稀疏的，所以图1的横坐标采用了对数坐标。表1中给出了打靶区间内的所有特征值，共有19个。Low和Lou在原始文章中给出了该序列的前两个值，他们给出的数值是0.425和2.55。
上述结果对应的参数设置如次：以1000的倍数作为分段点，将区间分为相应的10个
Figure 1. Variation of the error function in the shooting range of 10−1~104 for a2
图1. 误差函数在参数a2的打靶范围10−1~104内的变化曲线
Table 1. Eigenvalues defined by Low-Lou equation in the range of 0.1~10,000
表1. Low-Lou方程定义的特征值在0.1~10,000的数值结果
片段；每个片段均匀地插入3200个节点(含端点)，即，作为参数的“打靶点”；用龙格库塔方法求解初值问题时，特征函数的定义区间插入1999个节点。这些参数是通过试算确定的。下文结合物理背景验证和分析上述计算结果的正确性。
获得特征值后，可用龙格库塔方法计算对应的特征函数。特别地，Low-Lou给出了和的表达式：，。这意味着，只要得到特征值和特征函数，就可以按照公式(3)计算出空间任一点处的无力场。由于特征函数由数值方法得到，故而Low-Lou解法给出的无力场也称作“半解析解”。特别地，Wheatland等人 [9] 给出了衡量数值无力场的精确性的两个指标和，分别用于衡量数值磁场满足无力条件和无散条件的程度，其中“”表示所有数据取平均值。这两个指标是无量纲的标量，具体计算公式如下：
(9)
Table 2. Force-free and divergence-free merits for the Low-Lou solutions
表2. Low-Lou解的无力指标与无散指标
(10)
其中，表示节点附近小区域的体积，为该区域的表面积，表示节点处的电流，和分别表示节点处磁场和电流的大小。
注意，(9)中的本身是一种平均值，而(10)式中的是指节点处的局部值，由于数值磁场通常分布在立方网格的节点上，为此采用平均值来衡量数值磁场整体上满足无散条件的程度。本文采用Li等人 [10] 编写的程序计算这两个指标。在适当分辨率下，数值无力场的指标和应当分别达到（或小于）和的量级水平 [11] 。对应于表1所给的19个特征值，我们计算了相应的特征函数，并按公式(3)计算出19个Low-Lou解。计算中，我们统一采用立方网格，计算区域为。表2给出了这19个Low-Lou解的无力指标与无散指标。可以看到，两组指标都达到了满意的精度水平。特别地，对于序数较大的特征值，相应特征函数的零点增多，导致磁场结构更复杂，因而指标也略大。
4. 探讨与展望
非线性无力场是天体物理研究中的重要研究课题，而Low-Lou解法给出了求解轴对称非线性无力场的统一框架，它将问题转化为满足特定边界条件的带有未知参数的二阶非线性常微分方程。由于Low和Lou在原始文章中省略了技术细节，导致有关研究长期止步于最初求得的两组特解。本文提出的参数打靶法，能够有效地计算指定范围内的所有特征值，进而获得相应的数值无力场，可作为Low-Lou解法的补遗。需要指出的是，对于的情形，由于特征函数具有交叉型零点(zero-crossings)，方程(4)的非线性项会导致复数的出现，其物理意义不明显。为此，Flyer等人 [4] 考虑了方程(4)的变体形式，而Prasad等人 [7] 则引入了某种变换，兹不赘述。就的情形而言，可以验证，当特征值很大时，相应的数值无力场具有更复杂的几何与拓扑结构，预期可用于模拟太阳物理中的日珥和冕洞等重要物理现象。
参考文献