-
对许多数学和物理问题进行计算时，经常会得到无穷级数结果。对于一些收敛缓慢或发散的无穷级数，需要借助特殊的级数计算方法，才能得到有意义的数值结果。近几十年来，非线性数列变换的方法[1]逐渐发展，该方法能有效地加速收敛级数以及求发散级数的和，并在微扰论求解物理体系问题中得到越来越广泛的应用。
非谐振子是一种重要的物理系统，可以作为晶体振动的简化模型，且在量子场论中也有广泛的应用。在用微扰论求解强耦合条件下非谐振子基态能量本征值时，得到的解是发散级数，要计算无穷耦合极限必须借助特殊的求级数和技巧。非线性数列变换方法正适合解决这类问题。发散的级数经过非线性数列变换处理后，可以得到数学性质更好的新数列，从而收敛到一个广义的极限值（通常称为反极限），此过程即为发散级数求和。本文使用Weniger变换求非谐振子基态微扰级数和，从而计算四阶，六阶和八阶无穷耦合极限。
现有的求解非谐振子问题的工作多数采用浮点计算[2]，所得结果的精度都受到舍入误差的限制。本文利用计算机代数系统Maple的代数计算功能，进行完全有理化的计算，克服了舍入误差的限制。针对在计算机代数系统中进行高阶Weniger变换时遇到的内存迅速消耗问题，本文采用微扰级数系数分离计算的策略，并压缩非线性数列变换程序中关键数组的维数。通过这两种优化极大地解决了内存溢出的问题，得到精度极高的无穷耦合极限近似值。
-
谐振子是一种有精确解的物理模型，在许多物理问题中有广泛的应用。当它处于外场中时，就会产生附加的作用项，此时的体系称为非谐振子。非谐振子体系没有严格解，只能作近似处理，可以把非谐振子体系的哈密顿量写成[3]
其中β为耦合常数，
${\widehat{V}}^{(m)}$ 为非谐振子的新热能项${\widehat{H}}$ 。为谐振子的哈密顿量，为
式中
${\widehat{p}}^{\;2}$ 和${\widehat{x}}^{\;2}$ 分别为动能项和势能项，对于哈密顿量中的参数$ m $ 分别取2，3，4三个值，又可以将体系分为四阶，六阶和八阶非谐振子。按照Weniger的方法，本文将体系基态能量本征值的Schrödinger微扰级数形式解写成式（3）：
利用差分方程法[3]便可以逐项计算出该级数的所有系数
$ {b}_{n}^{\left(m\right)}\mathrm{。} $ Bender等[4]为了分析微扰级数发散速度构造了一种模型级数，对微扰级数的渐进性进行了详细的研究，结果表明m为2，3，4时非谐振子微扰级数的系数$ {b}_{n}^{\left(m\right)} $ 在$ n\to \mathrm{\infty } $ 处渐进展开式都是快速发散，从而微扰级数对应的模型级数会按照超几何级数发散，这些超几何级数对于所有耦合常数$ \beta $ 不为零的情况都是快速发散的。为了解决微扰级数发散速度过快的问题，对微扰级数进行重整化处理可以显著地降低发散速度，如Arteca等[5]对重整化各种算法做出了比较和总结。其中一种方式由Weniger得到[3]，他通过引入一对共轭变数：
得到了重整化后的哈密顿量
进一步根据变分法并结合其实际的物理意义，就可以得到一个表达式将待定常量
$ \tau $ 和耦合常数$ \beta $ 联系起来：
$ m= $ 2，3，4时分别对应$ {B}_{2}=3 $ ，$ {B}_{3}=45/4 $ ，$ {B}_{4}= 105/2 $ 。接下来引入重整化耦合常数
那么可以得到原哈密顿量中的原耦合常数
$ \beta $ 关于重整化后的耦合常数$ \kappa $ 的表达式
依照式（9）可以得到不同非谐振子体系中耦合常数变换的表达式。原耦合常数
$ \beta $ 符合物理意义的定义域是一个半无穷区间$ [0,\mathrm{\infty }) $ 。利用式（9），重整化可以将关于$ \beta $ 的半无穷区间投影到关于$ \kappa $ 的左闭右开单位区间$ \left[\mathrm{0,1}\right) $ 上。重整化方法将耦合常数的定义域重新构造成更小的区间，从而使发散微扰级数的求和变得相对简易，在耦合常数$ \beta $ 值较大时重整化的优势将尤为明显。本文将参量
$ {B}_{m} $ 的定义式（7）和新旧耦合常数间的关系式（9）带入重整化哈密顿量表达式中得到更简约的形式
其中
将式（10）带入Schrödinger方程则能得到重整化Schrödinger方程：
易看出重整化能量本征值
$ {E}_{R}^{\left(m\right)}\left(\kappa \right) $ 已由下标$ R $ 标记出来，此时耦合常数也变为$ \kappa $ 。基态下，非谐振子能量本征值
关于
$ \kappa $ 的级数展开为：
计算重整化级数系数
$ {c}_{n}^{\left(m\right)} $ 时本文采用和计算原级数系数$ {b}_{n}^{\left(m\right)} $ 类似的差分方程[3]来推导递推公式
其中
$ {G}_{j}^{\left(n\right)} $ ，$ 1\leqslant j\leqslant mn $ 均为计算级数系数$ {c}_{n}^{\left(m\right)} $ 的辅助参数。Banks等[6]将这两种哈密顿量解出的微扰级数系数作为特殊级数分析高阶渐进性，结果表明在指标
$ n\to \mathrm{\infty } $ 时级数系数$ {b}_{n}^{\left(m\right)} $ 是与一个表达式中含有$ {c}_{n}^{\left(m\right)} $ 的函数呈正相关，此外通过将重整化级数系数$ {c}_{n}^{\left(m\right)} $ 关于$ 1/n $ 做渐进展开并将首项和级数系数$ {b}_{n}^{\left(m\right)} $ 的首项作比较，最终还证明尽管重整化后的级数系数$ {c}_{n}^{\left(m\right)} $ 仍然关于$ n $ 呈阶乘增长趋势，但是增长速率相对于$ {b}_{n}^{\left(m\right)} $ 已经有了显著的减缓。依据哈密顿量（1）定义域的放大性质能得到，基态能
$ {E}^{\left(m\right)}\left(\beta \right) $ 的强耦合展开主项满足关系式[7]
其中
$ k_{m} $ 为无穷耦合极限。由式（12）和（15）可以看出重整化把基态能$ {E}^{\left(m\right)}\left(\beta \right) $ 变为重整化基态能$ {E}_{R}^{\left(m\right)}\left(\kappa \right) $ ，且$ \kappa =1 $ 对应$ \beta \to \mathrm{\infty } $ ，所以只需求无穷级数
就可以计算出无穷耦合极限
$k_{m} $ 。结合式（9）能得到$ k_{m} $ 的表达式
其中参数
$ {B}_{m} $ 的计算式已由式（7）给出。无穷耦合极限$ k_{m} $ 的计算所涉及的展开级数是重整化级数中最困难的问题，而有限场强$ \beta $ 对应耦合常数$ \kappa \in \left[\mathrm{0,1}\right) $ ，此时重整化微扰级数发散速度相对缓慢，求级数和更为容易。因此，计算无穷耦合极限的工作，对于其他任意有限场强中求解本征值问题有很好的借鉴意义。我们已知的最好的四阶，六阶和八阶非谐振子的无穷耦合极限结果是由Vinette[8]等通过重整化内投法（renormalized inner projection）计算到小数点后62,33，和21位小数。
-
非线性数列变换是一种加速级数收敛，或求发散级数和的方法。根据大量的实际计算结果[9-10]，Levin变换是效率最高且适用范围最广的一种非线性数列变换。然而，在求非谐振子基态本征值微扰级数和时，Levin变换所得结果出现了发散的现象，不能计算出高精度的非谐振子无穷耦合极限。考虑到Levin变换通常有很强的求发散级数和的能力，Weniger在保留Levin变换优点的基础上提出了Weniger变换[1]，该变换是以模型数列
为基础构造的，其中
$ {\omega }_{n} $ 为余项估计，待定参数$ \zeta =1 $ 。上式中采用了Pochhammer符号
如果（18）中的余项估计
$ {\omega }_{n} $ 与部分和之积能作为实际余项$ {r}_{n} $ 准确的近似式，那么Weniger变换就能给出高精度的极限$ s $ （发散级数的和通常称为反极限）近似值。利用差分算符的方法可以得到Weniger变换的定义式为
可以证明，可以采用三项递推公式
计算式（20）中的分子和分母。可以看出递推公式能在Weniger表中形成一个三角形。表中的元素
$ {J}_{k}^{\left(n\right)} $ 表示定义式（20）中分子或分母的计算值，如果用上标$ n $ 表示行标，下标$ k $ 表示列标，则可以构成一个呈左上三角形的二维矩阵：
其中变换阶
$ k $ 的最高阶$ {k}_{m}=l $ 。式（22）第一列中的元素是迭代计算过程的起始值，右边各列的元素都可以通过三项递推式（21）逐列计算得到。如果用
作为起始值，通过递推式（21）能得到Weniger变换式（20）的分子。如果用
作为起始值，则计算所得为Weniger变换式（20）的分母。将对应Weniger分子和分母表中相同位置的元素相除，则能得到Weniger变换
$ {\cal{I}}_{k}^{\left(n\right)} $ 结果表，显然这个表同式（22）一样是一个左上三角矩阵。根据分析[11]，四阶，六阶和八阶非谐振子重整化微扰级数属于各项正负交替的Stieltjes级数。因此，用余项估计
$ {\omega }_{n}=\mathrm{\Delta }{s}_{n}={a}_{n+1} $ 能有效地逼近实际余项，对应得到
在实际计算过程中，通常以待变数列部分和元素构成的有限子集
$ \{{s}_{1},{s}_{2},\cdots {s}_{l}\} $ 作为输入数据，我们称这个子集为待变序列。理论上，变换阶$ k $ 增加时产生的结果更精确。然而，在实际计算中，变换结果的精度开始时随着变换阶$ k $ 的增加而变高，在某一最佳变换阶得到有效位数最高的结果后，后续的结果就会逐渐变差。因此我们需要改变待变序列的长度$ l $ ，进行多次数列变换，进一步观察收敛的效果。本文借助计算机代数系统Maple实现了完全有理化的计算，数据在计算过程始终为有理数格式，所以不会产生舍入误差。然而，代价是单个数据的表示和计算会消耗更多的内存和机时。随着计算过程向高阶变换进行，Maple会消耗越来越多的内存，因此很容易发生内存溢出的问题。
针对内存溢出的问题，我们采用了两种优化方案。首先，根据（14）式可以看出，在计算微扰级数系数
$ {c}_{n}^{\left(m\right)} $ 时需要先计算辅助参数，且随着$ n $ 增大辅助参数的个数线性增加。将计算微扰级数系数过程和数列变换迭代过程程序分开，并把微扰级数系数输出到txt文件中，在计算变换结果前再以有理数格式读入，这样则能节省出可观的内存。其次，Weniger变换结果表是类似式（22）的三角形二维数组，如果把三项递推公式（21）计算的结果覆盖在原参数名上，则能将二维数组压缩成一维数组，节省出中间迭代过程消耗的内存。在迭代过程中，每列底部的数据会保持不变，因此我们最后所得一维数组正好是$ k=0,\mathrm{ }1,\mathrm{ }\cdots ,{k}_{m} $ 阶变换的结果，对应于式（22）矩阵中的反对角矩阵元。 -
使用上述的计算方法，我们用Weniger变换
$ t_{k}^{\left(n\right)}\left(\zeta ,{s}_{n}\right) $ 求出了四阶，六阶和八阶非谐振子基态无穷耦合极限。所有计算都是在Maple中完成。 -
在计算出微扰级数系数
$ {c}_{n}^{\left(m\right)} $ 后，由式（17）可以得到$ {E}_{R}^{\left(m\right)}\left(1\right) $ ，因此我们用部分和
作为输入数据。先以计算四阶非谐振子的无穷耦合极限
$k_{2} $ 为例，取$ m=2 $ 时四阶非谐振子微扰级数为$ {E}_{R}^{\left(2\right)}\left(1\right) $ ，表1列出了$ n=0\to 11 $ 的级数系数$ {c}_{n}^{\left(2\right)} $ 及其部分和数列$ {\sum }_{n}^{\left(2\right)} $ 的有理式结果。微扰级数系数$ {c}_{n}^{\left(2\right)} $ 及其部分和数列$ {\sum }_{n}^{\left(2\right)} $ 都是有理数格式参与运算，随着$ n $ 增大分子分母都会变长。表2给出了对应的浮点计算值，由此可以更直观地看出，微扰级数系数$ {c}_{n}^{\left(2\right)} $ 的绝对值随着$ n $ 增大而增大，级数和$ {E}_{R}^{\left(2\right)}\left(1\right) $ 是呈发散的趋势。$ n $ 
$ {c}_{n}^{\left(2\right)} $ 
$ {\sum }_{n}^{\left(2\right)} $ 
0 1 1 1 −1/4 3/4 2 -1/48 35/48 3 1/64 143/192 4 −791/27648 19801/27648 5 7273/110592 86477/110592 6 −1462711/7962624 4763633/7962624 7 19238731/31850496 1418269/1179648 8 −20961986815/9172942848 −9933527071/9172942848 9 119587982023/12230590464 319029837785/36691771392 10 −123340039323181/
2641807540224−100369891002661/
264180754022411 2600833240032557/
105672301608962199353676021913/
10567230160896表 1 四阶非谐振子微扰级数系数
$ {c}_{n}^{\left(2\right)} $ 及其部分和数列$ {\sum }_{n}^{\left(2\right)} $ Table 1. Coefficients
$ {c}_{n}^{\left(2\right)} $ and partial sums$ {\sum }_{n}^{\left(2\right)} $ of perturbation series for quartic anharmonic oscillator$ n $ 
$ {c}_{n}^{\left(2\right)} $ 
$ {\sum }_{n}^{\left(2\right)} $ 
0 1.0000000000 1.0000000000 1 −0.2500000000 0.7500000000 2 −0.0208333333 0.7291666667 3 0.0156250000 0.7447916667 4 −0.0286096644 0.7161820023 5 0.0657642506 0.7819462529 6 −0.1836971079 0.5982491451 7 0.6040323830 1.2022815280 8 −2.2851975819 −1.0829160538 9 9.7777766638 8.6948606099 10 −46.6877459638 −37.9928853539 11 246.1225127524 208.1296273985 表 2 null的浮点数计算结果
Table 2. Floating-type coefficients
$ {c}_{n}^{\left(2\right)} $ and partial sums$ {\sum }_{n}^{\left(2\right)} $ ofperturbation series for quartic anharmonic oscillator我们在实际计算时总是从
$ {\sum }_{1}^{\left(m\right)} $ 开始,用形如
的待变序列做非线性数列变换。四阶非谐振子
$ m=2 $ ，若取$ l=10 $ 得到待变序列$ \{{\sum }_{1}^{\left(2\right)},{\sum }_{2}^{\left(2\right)},\cdots {\sum }_{10}^{\left(2\right)}\} $ ，对应着表1中第三列$ n=1\to 10 $ 的10个待变序列元素。对$ \{{\sum }_{1}^{\left(2\right)},{\sum }_{2}^{\left(2\right)},\cdots {\sum }_{10}^{\left(2\right)}\} $ 做Weniger变换$ t_{l}^{\left(n\right)}\left(\zeta ,{s}_{n}\right) $ ，即公式（25），得到级数和$ {E}_{R}^{\left(2\right)}\left(1\right) $ 的各阶数列变换结果为$ t_{k}^{\left(10-k\right)} $ ，$ k=0\to 9 $ 。根据式（17），将变换后所得级数和乘一个常数则能算出各阶非谐振子的无穷耦合极限$ k_{m} $ ，四阶非谐振子无穷耦合极限近似值为
计算结果列于表3第二列。由表3可以看出，无穷耦合极限
$ k_{2} $ 的吻合的有效位数开始时随着变换阶$ k $ 的增加而变多，到了第7阶后有效位数又开始减少。$ k $ 
$ k_{2} $ 
$ {\delta }_{k} $ 
0 −5.47952225763255·101 1 1.80372780913117 0.129 2 1.06589749248284 0.132 3 1.06038019425837 2.258 4 1.06035927005928 4.679 5 1.06036167836480 5.618 6 1.06036202937495 6.455 7 1.06036191647094 6.947 8 1.06036233326828 6.380 9 1.06036453736592 5.657 表 3 序列（27）使用Weniger变换的结果，第二列是变换后由（28）得到的四阶无穷耦合极限
$k_{2}$ 近似值Table 3. Numerical results of Weniger transformation for the string (27), and entries in the second column indicates approximations of quartic infinite coupling limit
$k_{2}$ according to (28)为了更定量且直观地考察Weniger变换过程的收敛情况，定义：
$ {\delta }_{k} $ 表示相邻两次迭代结果的符合程度，大致对应于两个十进制浮点数取得一致的位数。不同变换阶$ k $ 时计算结果的$ {\delta }_{k} $ 不同，通常$ {\delta }_{k} $ 值越大对应的近似值精度也就越高。从表3可以看出，通过数列变换的方法，我们从发散级数所包含的数学信息中提取出了收敛的有限值，原本发散的微扰级数因此获得了物理意义。$ {\delta }_{k} $ 先随着变换阶$ k $ 的增加而不断增加，在最佳变换阶$ k=7 $ 时最大，然后转为变小。因此最高精度为$ {\Delta }_{10}=6.947 $ ，对应收敛到7位有效数字的结果为$ k_{2}=1.060\;361 $ 。为了得到精度更高的计算结果，必须使用更长的待变序列（27）进行更高阶的Weniger变换。我们选取序列
$ \{{\sum }_{1}^{\left(2\right)},{\sum }_{2}^{\left(2\right)},\cdots {\sum }_{100}^{\left(2\right)}\} $ ，以这个序列作为输入数据可以计算出最高阶$ {k}_{m}=99 $ 的Weniger变换结果，如表4所示。可以看到$ k=88 $ 阶变换时，最大的$ {\delta }_{k}=29.000 $ ，对应的最佳收敛值$ k_{2}=1.060\;362\;090\;484\; 182\;899\;647\;046\;016\;70 $ ，有效数字可取到小数点后大致第29位。为了比较取不同长度$ l $ 待变序列（27）的变换结果（即有效数字的位数），我们把这时的$ {\delta }_{k} $ 值同时定义为$ {\Delta }_{l} $ 的值：$ {\Delta }_{100}={\delta }_{88}=29.000 $ 。$ k $ 
$ k_{2} $ 
$ {\delta }_{k} $ 
83 1.060 362 090 484 182 899 647 046 007 50 24.028 84 1.060 362 090 484 182 899 647 046 024 36 25.773 85 1.060 362 090 484 182 899 647 046 019 00 26.271 86 1.060 362 090 484 182 899 647 046 016 95 26.688 87 1.060 362 090 484 182 899 647 046 016 69 27.585 88 1.060 362 090 484 182 899 647 046 016 70 29.000 89 1.060 362 090 484 182 899 647 046 016 34 27.444 90 1.060 362 090 484 182 899 647 046 014 87 26.833 91 1.060 362 090 484 182 899 647 046 011 57 26.481 92 1.060 362 090 484 182 899 647 046 008 74 26.548 93 1.060 362 090 484 182 899 647 046 019 40 25.972 表 4
$ l=100 $ 时Weniger变换的结果，$ k=88 $ 时$ {\delta }_{k}=29.000 $ Table 4. Numerical results of Weniger transformation for
$ l=100 $ , and$ {\delta }_{k}=29.000 $ when$ k=88 $ 为了更深入地考查Weniger变换求四阶非谐振子微扰级数和的能力，我们按照与
$ l=100 $ 时相同的方法，依次计算了待变序列长度$ l=200,\mathrm{ }300,\mathrm{ }\dots ,\mathrm{ }800 $ 的Weniger变换结果，整理成表5。计算Weniger变换所选择的待变序列越长，所需要的微扰级数系数$ {c}_{n}^{\left(2\right)} $ 越多，可以进行更高阶的变换，给出的结果精度越高。最高精度$ {\Delta }_{l}=105.646 $ ，是由待变序列长度$ l=800 $ 时，变换到$ k=756 $ 得到l Δl l Δl 100 29.000 500 77.658 200 43.602 600 89.463 300 56.339 700 96.025 400 68.265 800 105.646 表 5 对不同长度l待变序列（27）进行Weniger变换得到无穷耦合极限
$ k_{2} $ 的有效位数$ {\Delta }_{l} $ Table 5. The number of decimal significant digits
$ {\Delta }_{l} $ of approximationsof infinite coupling limit$ k_{2} $ given by Weniger transformationfor various lengths of strings (27)
比现有的最准确的内投法结果[8]还要高出大约43位有效数字。
我们使用图形软件Origin，以拟合函数
$ {\Delta }_{l}=a+bl+c{l}^{2} $ 对表5中的$ {\Delta }_{l} $ 和$ l $ 数据做二次多项式拟合，可以得到图1。其中，截距$ a=14.561\;07\pm 1.315\;21 $ ，一阶系数$ b=0.153\;49\pm \mathrm{ }0.006\;71 $ ，二阶系数$ \mathrm{c}=-5.0519\times {10^{ - 5}}\pm 7.27312\times {10^{ - 6}} $ ，残差平方和为$ 4.443\;45 $ ，R2(COD)为$ 0.999\;1 $ 。从图1可以看出吻合的有效数字位数$ {\Delta }_{l} $ 关于待变序列长度$ l $ 大致是线性增加的，这说明我们只需增加待变序列长度就可以计算出有效数字更多的四阶非谐振子耦合极限$ k_{2} $ 。
图 1 四阶非谐振子无穷耦合极限
$ k_{2} $ 近似值的有效位数$ {\Delta}_{l} $ 与待变序列式（27）长度$ l $ 之间的关系。有效位数$ {\Delta}_{l} $ 是由Weniger变换$ t_{l}^{\left(1\right)}\left(1,{\sum }_{1}^{\left(2\right)}\right) $ 结果计算所得。拟合函数为$ {\Delta}_{l}=a+bl+c{l}^{2} $ Figure 1. Dependence of the number of decimal significant digits
$ {\Delta}_{l} $ on lengths$ l $ of strings to be transformed formula (27).$ {\Delta}_{l} $ is given by Weniger transformation$ t_{l}^{\left(1\right)}\left(1,{\sum }_{1}^{\left(2\right)}\right) $ . Fitting function is$ {\Delta}_{l}=a+bl+c{l}^{2} $ -
六阶非谐振子微扰级数系数增长更迅速，文献[3]中在指标
$ n $ 相同时微扰系数$ {c}_{n}^{\left(3\right)} $ 相对于四阶非谐振子的系数$ {c}_{n}^{\left(2\right)} $ 更大，微扰级数发散更快。因此待变序列式（27）的长度相同时，六阶非谐振子的计算量更大。依照四阶非谐振子的情形，我们使用Weniger变换求微扰级数和$ {E}_{R}^{\left(3\right)}\left(1\right) $ ，由式（17）也可以求出六阶非谐振子无穷耦合极限$ k_{3} $ 近似值。以$\{{\sum }_{1}^{\left(3\right)},{\sum }_{2}^{\left(3\right)},\cdots {\sum }_{100}^{\left(3\right)}\}$ 做第一次变换，初步观察收敛效果，变换结果见表6。可以看出，在变换阶$ k=92 $ 时，六阶非谐振无穷子耦合极限$k_{3} $ 的有效位数$ {\Delta}_{100}=11.053 $ 。$ k $ 
$ k_{3} $ 
$ {\delta }_{k} $ 
87 1.144 802 440 193 7.866 88 1.144 802 446 940 8.171 89 1.144 802 453 804 8.163 90 1.144 802 453 891 10.063 91 1.144 802 453 780 9.955 92 1.144 802 453 788 11.053 93 1.144 802 453 801 10.896 94 1.144 802 453 816 10.815 95 1.144 802 453 850 10.471 96 1.144 802 453 903 10.279 97 1.144 802 453 954 10.295 表 6
$ l=100 $ 时Weniger变换的结果，$ k=92 $ 时$ {\delta }_{k}=11.053 $ Table 6. Numerical results of Weniger transformation for
$ l= 100 $ , and$ {\delta }_{k}=11.053 $ when$ k=92 $ 我们依次对长度
$ l=200,\mathrm{ }300,\mathrm{ }\dots ,\mathrm{ }1\;000 $ 的待变序列进行了Weniger变换。我们服务器内存最高可以负载待变序列长度$ l=1\;000 $ 时的计算量，具体数据列于表7。表7显示，变换$ t_{1000}^{\left(1\right)}\left(1,{\sum }_{1}^{\left(3\right)}\right) $ 的最高精度为$ {\Delta}_{l}=32.711 $ ，即计算结果收敛到小数点后33位十进制数字$ k $ 
$ k_{3} $ 
$ {\delta }_{k} $ 
963 1.144 802 453 797 052 763 765 457 535 812 934 26.779 964 1.144 802 453 797 052 763 765 457 534 188 410 26.789 965 1.144 802 453 797 052 763 765 457 534 148 665 28.401 966 1.144 802 453 797 052 763 765 457 534 149 462 30.099 967 1.144 802 453 797 052 763 765 457 534 149 546 31.076 968 1.144 802 453 797 052 763 765 457 534 149 544 32.711 969 1.144 802 453 797 052 763 765 457 534 149 536 32.102 970 1.144 802 453 797 052 763 765 457 534 149 524 31.931 971 1.144 802 453 797 052 763 765 457 534 149 513 31.973 972 1.144 802 453 797 052 763 765 457 534 149 519 32.230 表 7
$ l=1\;000 $ 时Weniger变换的结果，$ k=968 $ 时$ {\delta }_{k}=32.711 $ Table 7. Numerical results of Weniger transformation for
$ l=1\;000 $ , and$ {\delta }_{k}=32.711 $ when$ k=968 $ 
和内投法结果[8]
基本吻合，只是最后一位带误差的小数不相同。如果按照图1的方法绘制出拟合曲线，也能得到一个大致的线性关系。如能继续增加待变序列长度，则能得到比内投法更多的有效位数。
-
以长度
$ l=100 $ 待变序列（27）经过Weniger变换所得八阶非谐振子无穷耦合极限$ k_{4} $ 结果如表8所示。可以看出，四阶，六阶和八阶非谐振子在待变序列长度$ l=100 $ 时，有效位数$ {∆}_{l} $ 逐渐减少，依次为29，11和6。对于同样长度的待变序列，由于级数的发散特性更加明显，$ m $ 值大阶数高的非谐振子需要消耗更大的计算量和内存空间，却只能得到更差的结果。$ k $ 
$ k_{4} $ 
$ {\delta }_{k} $ 
92 1.225 754 3 4.081 93 1.225 810 1 4.995 94 1.225 821 4 5.515 95 1.225 818 9 5.523 96 1.225 815 9 5.552 97 1.225 813 1 5.785 98 1.225 811 4 5.962 99 1.225 812 3 5.283 表 8
$ l=100 $ 时Weniger变换的结果，$ k=98 $ 时$ {\delta }_{k}=5.962 $ Table 8. Numerical results of Weniger transformation for
$ l=100 $ , and$ {\delta }_{k}=5.962 $ when$ k=98 $ 我们依次对长度
$ l=200,\mathrm{ }300,\mathrm{ }\dots ,\mathrm{ }900 $ 的待变序列（27）进行了Weniger变换。在待变序列长度$ l>500 $ 后，出现了服务器内存溢出的问题。我们通过优化数组结构等方法，成功计算出了待变序列长度$ l=900 $ 时八阶非谐振子的无穷耦合极限，结果见表9。从表9可以看出变换阶$ k=890 $ 时得到最高有效位数$ {\Delta}_{900}= 13.488 $ ，即第890阶变换结果和第889阶变换结果大致有14位吻合的有效数字，这时已到达我们计算资源的极限。$ k $ 
$ k_{4} $ 
$ {\delta }_{k} $ 
887 1.225 820 113 800 842 11.287 888 1.225 820 113 800 340 12.300 889 1.225 820 113 800 514 12.759 890 1.225 820 113 800 547 13.488 891 1.225 820 113 800 485 13.209 892 1.225 820 113 800 328 12.803 893 1.225 820 113 800 002 12.488 894 1.225 820 113 799 447 12.255 895 1.225 820 113 798 641 12.094 896 1.225 820 113 797 630 11.996 表 9
$ l=900 $ 时Weniger变换的结果，$ k=890 $ 时$ {\delta }_{k}=13.488 $ Table 9. Numerical results of Weniger transformation for
$ l=900 $ , and$ {\delta }_{k}=13.488 $ when$ k=890 $ -
本文研究了Weniger变换求四阶，六阶和八阶非谐振子无穷耦合极限的应用。计算结果表明，Weniger变换有很好的求发散级数和能力。在计算四阶非谐振子耦合极限
$ k_{2} $ 时，增加待变序列$\{{\sum }_{1}^{\left(2\right)}, {\sum }_{2}^{\left(2\right)},\cdots {\sum }_{l}^{\left(2\right)}\}$ 长度$ l $ ，进行更高阶变换，总能得到迅速收敛的结果。在待变序列长度$ l=800 $ 时，最高精度$ {\Delta}_{l}=105.646 $ ，这个结果比我们已知的通过内投法得到的最准确结果还要高出大约43位有效数字。我们对无穷耦合极限的计算结果与待变序列的长度（$ l $ 值）之间的关系进行了系统的研究，得出了图1所示的结果：有效位数$ {\Delta}_{l} $ ~$ l $ 大致是一个线性关系。由此可以看出，只要有足够的计算资源，我们就可以通过增加待变序列的长度，进行相应的非线性变换就能得到更精确的，有效位数大致线性增长的计算结果。在Maple计算过程中，每个系数是巨大的有理分数形式，不仅占据了大量内存，还会明显增加计算量。六阶和八阶非谐振子微扰级数系数比四阶的增长更快，在进行Weniger变换时需要更多的内存和计算量，计算无穷耦合极限也就更为困难。我们计算无穷耦合极限
$k_{3} $ 和$k_{4} $ 时遇到了服务器内存溢出的问题，为此我们改进了Weniger的程序数组结构，将二维数组通过覆盖的编程技巧将其压缩为一维数组。一维数组的优化方案节省了巨大的内存，这个改进使我们可以计算出比Weniger更准确的结果。不然的话，即使是在我们现有的计算条件下也很难得出比Weniger [3]更好的计算结果。此外，由公式（14）可以看出，将计算微扰级数系数$ {c}_{n}^{\left(m\right)} $ 部分的程序分离出迭代的变换过程，也能节省出大量保存辅助参数时被浪费的内存。微扰级数的微扰阶越高，待变序列长度越长，这种优化方案的优势就更显著。表10中汇总了我们优化方案的计算结果，以及和其他方法所得结果的比较。表 10 四阶，六阶和八阶无穷耦合极限计算结果与其他方法的比较
Table 10. Comparison between present numerical results of quartic,sextic and octic infinite coupling limits and other methods
如果有更多的计算机内存资源，则可以计算出更多位数的无穷耦合极限。而且当内存足够充足时，若使用Maple的并行运算功能，将极大地提升计算速度，节省CPU时间。因此，这里内存资源则是更为关键因素。
微扰论是求解许多物理体系的常用近似计算方法。实际计算得到的无穷级数为缓慢收敛和发散的情况也是非常普遍的。于是，通过非线性数列变换改善无穷级数的收敛性质具有普遍的实际意义[12]。
Weniger变换在求非谐振子无穷耦合极限中的应用
Application of the Weniger’s Transformation in Approximations of Infinite Coupling Limits of Anharmonic Oscillators
-
摘要: 研究了Weniger变换求非谐振子基态能强耦合微扰展开发散级数和，并计算出无穷耦合极限。使用计算机代数系统Maple克服了舍入误差对数值计算的负面影响，代价是每个数据的表示和运算会消耗更多的内存。提出一种优化数组结构的方案，有效地缓解了内存压力，在现有的内存资源下得到高精度的计算结果。Abstract: During decades, nonlinear sequence transformations method has been well developed in fields of mathematics and physics, and extensive simulation results have demonstrated its power of the acceleration of convergence and the summation of divergent series. The perturbation expansions for the infinite coupling limits of the quartic, sextic and octic anharmonic oscillators are strongly divergent, and renormalization techniques shall be used to slow down its rate of divergence. This paper presents the performance of Weniger’s transformation in summation of the renormalized perturbation series, and gives numerical results of infinite coupling limits. With the help of computer algebra system Maple, which has abilities of rational arithmetics, we can get rid of the bad effect of rounding errors. However, Maple consumes large amounts of memory resources to store data and calculate, as a result memory overflow occurs frequently. Aiming at the above problem, this paper proposes a method to compress the dimensions of arrays in order to reduce load of storage, and thus we can obtain more accurate approximations of infinite coupling limits than the known method.
-
图 1 四阶非谐振子无穷耦合极限
$ k_{2} $ 近似值的有效位数$ {\Delta}_{l} $ 与待变序列式（27）长度$ l $ 之间的关系。有效位数$ {\Delta}_{l} $ 是由Weniger变换$ t_{l}^{\left(1\right)}\left(1,{\sum }_{1}^{\left(2\right)}\right) $ 结果计算所得。拟合函数为$ {\Delta}_{l}=a+bl+c{l}^{2} $ Figure 1. Dependence of the number of decimal significant digits
$ {\Delta}_{l} $ on lengths$ l $ of strings to be transformed formula (27).$ {\Delta}_{l} $ is given by Weniger transformation$ t_{l}^{\left(1\right)}\left(1,{\sum }_{1}^{\left(2\right)}\right) $ . Fitting function is$ {\Delta}_{l}=a+bl+c{l}^{2} $ 表 1 四阶非谐振子微扰级数系数
$ {c}_{n}^{\left(2\right)} $ 及其部分和数列$ {\sum }_{n}^{\left(2\right)} $ Table 1. Coefficients
$ {c}_{n}^{\left(2\right)} $ and partial sums$ {\sum }_{n}^{\left(2\right)} $ of perturbation series for quartic anharmonic oscillator$ n $ $ {c}_{n}^{\left(2\right)} $ $ {\sum }_{n}^{\left(2\right)} $ 0 1 1 1 −1/4 3/4 2 -1/48 35/48 3 1/64 143/192 4 −791/27648 19801/27648 5 7273/110592 86477/110592 6 −1462711/7962624 4763633/7962624 7 19238731/31850496 1418269/1179648 8 −20961986815/9172942848 −9933527071/9172942848 9 119587982023/12230590464 319029837785/36691771392 10 −123340039323181/
2641807540224−100369891002661/
264180754022411 2600833240032557/
105672301608962199353676021913/
10567230160896
下载: 导出CSV
表 2 null的浮点数计算结果
Table 2. Floating-type coefficients
$ {c}_{n}^{\left(2\right)} $ and partial sums$ {\sum }_{n}^{\left(2\right)} $ ofperturbation series for quartic anharmonic oscillator$ n $ $ {c}_{n}^{\left(2\right)} $ $ {\sum }_{n}^{\left(2\right)} $ 0 1.0000000000 1.0000000000 1 −0.2500000000 0.7500000000 2 −0.0208333333 0.7291666667 3 0.0156250000 0.7447916667 4 −0.0286096644 0.7161820023 5 0.0657642506 0.7819462529 6 −0.1836971079 0.5982491451 7 0.6040323830 1.2022815280 8 −2.2851975819 −1.0829160538 9 9.7777766638 8.6948606099 10 −46.6877459638 −37.9928853539 11 246.1225127524 208.1296273985
下载: 导出CSV
表 3 序列（27）使用Weniger变换的结果，第二列是变换后由（28）得到的四阶无穷耦合极限
$k_{2}$ 近似值Table 3. Numerical results of Weniger transformation for the string (27), and entries in the second column indicates approximations of quartic infinite coupling limit
$k_{2}$ according to (28)$ k $ $ k_{2} $ $ {\delta }_{k} $ 0 −5.47952225763255·101 1 1.80372780913117 0.129 2 1.06589749248284 0.132 3 1.06038019425837 2.258 4 1.06035927005928 4.679 5 1.06036167836480 5.618 6 1.06036202937495 6.455 7 1.06036191647094 6.947 8 1.06036233326828 6.380 9 1.06036453736592 5.657
下载: 导出CSV
表 4
$ l=100 $ 时Weniger变换的结果，$ k=88 $ 时$ {\delta }_{k}=29.000 $ Table 4. Numerical results of Weniger transformation for
$ l=100 $ , and$ {\delta }_{k}=29.000 $ when$ k=88 $ $ k $ $ k_{2} $ $ {\delta }_{k} $ 83 1.060 362 090 484 182 899 647 046 007 50 24.028 84 1.060 362 090 484 182 899 647 046 024 36 25.773 85 1.060 362 090 484 182 899 647 046 019 00 26.271 86 1.060 362 090 484 182 899 647 046 016 95 26.688 87 1.060 362 090 484 182 899 647 046 016 69 27.585 88 1.060 362 090 484 182 899 647 046 016 70 29.000 89 1.060 362 090 484 182 899 647 046 016 34 27.444 90 1.060 362 090 484 182 899 647 046 014 87 26.833 91 1.060 362 090 484 182 899 647 046 011 57 26.481 92 1.060 362 090 484 182 899 647 046 008 74 26.548 93 1.060 362 090 484 182 899 647 046 019 40 25.972
下载: 导出CSV
表 5 对不同长度l待变序列（27）进行Weniger变换得到无穷耦合极限
$ k_{2} $ 的有效位数$ {\Delta }_{l} $ Table 5. The number of decimal significant digits
$ {\Delta }_{l} $ of approximationsof infinite coupling limit$ k_{2} $ given by Weniger transformationfor various lengths of strings (27)l Δl l Δl 100 29.000 500 77.658 200 43.602 600 89.463 300 56.339 700 96.025 400 68.265 800 105.646
下载: 导出CSV
表 6
$ l=100 $ 时Weniger变换的结果，$ k=92 $ 时$ {\delta }_{k}=11.053 $ Table 6. Numerical results of Weniger transformation for
$ l= 100 $ , and$ {\delta }_{k}=11.053 $ when$ k=92 $ $ k $ $ k_{3} $ $ {\delta }_{k} $ 87 1.144 802 440 193 7.866 88 1.144 802 446 940 8.171 89 1.144 802 453 804 8.163 90 1.144 802 453 891 10.063 91 1.144 802 453 780 9.955 92 1.144 802 453 788 11.053 93 1.144 802 453 801 10.896 94 1.144 802 453 816 10.815 95 1.144 802 453 850 10.471 96 1.144 802 453 903 10.279 97 1.144 802 453 954 10.295
下载: 导出CSV
表 7
$ l=1\;000 $ 时Weniger变换的结果，$ k=968 $ 时$ {\delta }_{k}=32.711 $ Table 7. Numerical results of Weniger transformation for
$ l=1\;000 $ , and$ {\delta }_{k}=32.711 $ when$ k=968 $ $ k $ $ k_{3} $ $ {\delta }_{k} $ 963 1.144 802 453 797 052 763 765 457 535 812 934 26.779 964 1.144 802 453 797 052 763 765 457 534 188 410 26.789 965 1.144 802 453 797 052 763 765 457 534 148 665 28.401 966 1.144 802 453 797 052 763 765 457 534 149 462 30.099 967 1.144 802 453 797 052 763 765 457 534 149 546 31.076 968 1.144 802 453 797 052 763 765 457 534 149 544 32.711 969 1.144 802 453 797 052 763 765 457 534 149 536 32.102 970 1.144 802 453 797 052 763 765 457 534 149 524 31.931 971 1.144 802 453 797 052 763 765 457 534 149 513 31.973 972 1.144 802 453 797 052 763 765 457 534 149 519 32.230
下载: 导出CSV
表 8
$ l=100 $ 时Weniger变换的结果，$ k=98 $ 时$ {\delta }_{k}=5.962 $ Table 8. Numerical results of Weniger transformation for
$ l=100 $ , and$ {\delta }_{k}=5.962 $ when$ k=98 $ $ k $ $ k_{4} $ $ {\delta }_{k} $ 92 1.225 754 3 4.081 93 1.225 810 1 4.995 94 1.225 821 4 5.515 95 1.225 818 9 5.523 96 1.225 815 9 5.552 97 1.225 813 1 5.785 98 1.225 811 4 5.962 99 1.225 812 3 5.283
下载: 导出CSV
表 9
$ l=900 $ 时Weniger变换的结果，$ k=890 $ 时$ {\delta }_{k}=13.488 $ Table 9. Numerical results of Weniger transformation for
$ l=900 $ , and$ {\delta }_{k}=13.488 $ when$ k=890 $ $ k $ $ k_{4} $ $ {\delta }_{k} $ 887 1.225 820 113 800 842 11.287 888 1.225 820 113 800 340 12.300 889 1.225 820 113 800 514 12.759 890 1.225 820 113 800 547 13.488 891 1.225 820 113 800 485 13.209 892 1.225 820 113 800 328 12.803 893 1.225 820 113 800 002 12.488 894 1.225 820 113 799 447 12.255 895 1.225 820 113 798 641 12.094 896 1.225 820 113 797 630 11.996
下载: 导出CSV
-
[1] WENIGER E J. Nonlinear sequence transformations for the acceleration of convergence and the summation of divergent series[J]. Computer Physics Reports, 1989, 10(5): 189-371. [2] KILLINGBECK J. The harmonic oscillator with λx.M perturbation[J]. Journal of Physics A: Mathematical and Theoretical, 1980, 13(1): 49-56. [3] WENIGER E J, ČÍŽEK J, VINETTE F. The summation of the ordinary and renormalized perturbation series for the ground state energy of the quartic, sextic, and octic anharmonic oscillators using nonlinear sequence transformations[J]. Journal of Mathematical Physics, 1992, 34(2): 571-609. [4] BENDER C M, WU T T. Large-order behavior of perturbation theory[J]. Physics Review Letters, 1971, 27(7): 461-465. [5] ARTECA G A, FERNÁNDEZ F M, CASTRO E A. Large Order Perturbation Theory and Summation Methods in Quantum Mechanics[M]. Berlin: Springer, 1990. [6] BANKS T I, BENDER C M. Anharmonic oscillator with polynomial self-interaction[J]. Journal of Mathematical Physics, 1972, 13(9): 1320-1324. [7] SIMON B. In Cargèse Lectures in Physics, edited by BESSIS D[M]. New York: Gordon and Breach, 1972. [8] VINETTE F, ČÍŽEK J. Upper and lower bounds of the ground state energy of anharmonic oscillators using renormalized inner projection[J]. Journal of Mathematical Physics, 1991, 32(12): 3392-3404. [9] SMITH D A, FORD W F. Acceleration of linear and logarithmic convergence[J]. Society for Industrial Applied Mathematics Journal on Numerical Analysis, 1979, 16(2): 223-240. [10] SMITH D A, FORD W F. Numerical comparisons of nonlinear convergence accelerators[J]. Mathematics of Computation, 1982, 38(158): 481-499. [11] SIMON B. Coupling constant analyticity for the anharmonic oscillator[J]. Annals of Physics, 1970, 58(1): 76-136. [12] BELKIČ D. New hybrid non-linear transformations of divergent perturbation series for quadratic Zeeman effects[J]. Journal of Physics A: Mathematical and Theoretical, 1989, 22(15): 3003-3010. -
点击查看大图
计量
- 文章访问数: 316
- HTML全文浏览量: 200
- PDF下载量: 1
- 被引次数: 0