一个代数方程根的级数解问题

数学星空

求下列关于\(x\)方程根的近似解：

1. \(\frac{1}{2x}+\frac{1}{x-1}+\frac{1}{x-2}+\dots+\frac{1}{x-n}=a\)

2. \(\frac{1}{x}+\frac{2x}{x^2-1}+\frac{2x}{x^2-2^2}+\dots+\frac{2x}{x^2-n^2}=a\)


根据相关资料已有结果(对于最大根\(x_n\)的求解)：

对1. \(x\approx\frac{n+\frac{1}{2}}{1-e^{-a}}\)

对2. \(x\approx\frac{(n+\frac{1}{2})(1+e^{-a})}{1-e^{-a}}\)

我们现在需要做的是给出更精确的级数解。

282842712474

对于1，设$x > n$，则可以用欧拉－麦克劳林求和公式(https://en.wikipedia.org/wiki/Euler%E2%80%93Maclaurin_formula)
$$a=-\frac{1}{2x}+\sum_{k=0}^n \frac{1}{x-k}=-\frac{1}{2x}+\int_0^n \frac{1}{x-k}dk+\frac{1}{2}\left(\frac{1}{x-1}+\frac{1}{x-n}\right)+\dots=-\frac{1}{2x}+\ln x-\ln(x-n)+\frac{1}{2}\left(\frac{1}{x-1}+\frac{1}{x-n}\right)+\dots$$
但是，如果考虑级数展开，需要分情况，$a$很大的时候，和$a$很小的时候。

282842712474

考虑$a$很小时，也就是$x$足够大时，那么就有近似
$$a=\ln\left(\frac{x}{x-n}\right)$$
得到
$$x=\frac{n}{1-e^{-a}}$$
略有不一样

kastin

这里有个逆级数链接 https://en.wikipedia.org/wiki/Lagrange_inversion_theorem，不过由于函数是分式，若用幂级数的话，`n`和`a`的取值对其导函数影响较大，有点困难。

1. 先从形式上来看，函数在整点处存在极点，有点类似余切函数，我们来尝试下

1. f[x_, n_] := 1/2/x + Sum[1/(x - k), {k, 1, n}];
   
2. sol = NSolve[D[f[x, 6], {x, 2}] == 0, x, Reals];
   
3. data = {x, f[x, 6]} /. sol (* 拐点坐标 *)
   
4. Plot[{f[x, 6],  Pi Cot[Pi x]}, {x, 0, 10}, PlotStyle -> {Blue, Gray},
   
5. Epilog -> {Orange, Line[data], PointSize[Medium], Red, Point[data]},
   
6. PlotLegends -> "Expressions"]

复制代码

问题1

图中红色的点为拐点，上与余切函数的拐点横坐标很接近：
{{x -> 0.43739}, {x -> 1.49746}, {x -> 2.49943}, {x -> 3.50012}, {x -> 4.50095}, {x -> 5.50411}}
很明显，图像关于 `x=\lceil n/2 \rceil` 对称。
若 `a` 的绝对值取值越大（比如大于10），那么求得的根可以用反余切函数来近似。共有 `n+1` 个根，当然精确的根分布并不具有周期性。近似根越是靠近对称轴的则越精确。
令 `x_k\;(k=1,2,\cdots,n+1)` 为原方程从小到大排列的第 `k` 个根的近似值，那么有$$\pi \cot \pi x= a\;(|a|>10)$$因此$$x_k=k+\frac{1}{\pi}\mathrm{arccot}\frac{a}{\pi}$$

2. 仿照1的方法，仍然还是余切函数逼近，下面是图像

1. g[x_, n_] := 1/x + Sum[2 x/(x^2 - k^2), {k, 1, n}];
   
2. sol = NSolve[D[g[x, 6], {x, 2}] == 0, x, Reals]
   
3. data = {x, g[x, 6]} /. sol;
   
4. Plot[{g[x, 6],  Pi  Cot[Pi x]}, {x, -7, 7}, PlotStyle -> {Blue, Gray},
   
5.   Epilog -> {Orange, Line[data], PointSize[Medium], Red,
   
6.    Point[data]}, PlotLegends -> "Expressions"]

复制代码

问题2

红色的点为拐点，同样与余切拐点横坐标很接近。
{{x->-5.50425},{x->-4.50117},{x->-3.5005},{x->-2.50025},{x->-1.50012},{x->-0.500037},{x->0.500037},{x->1.50012},{x->2.50025},{x->3.5005},{x->4.50117},{x->5.50425}}
可以看出，问题2的图像是关于y轴(`x=0`)对称分布的，共有 `2n+1` 个根。若 `a` 的绝对值取值越大（比如大于10），那么求得的根可以用反余切函数来近似。近似根越是靠近对称轴（`x=0`）则越精确。
求解方法与上面雷同，就不赘述了。

从上面可以看出，如果 `|a|<10`，越是远离对称轴，近似根精确性越差，但是目前没想到什么好方法。

补充内容 (2016-7-6 12:55):
更正一下，应该是“近似反对称”。仔细看的话，图像并不是完全反对称的。上面的近似解与 `n` 的取值影响不大，反倒是 `a` 的取值影响更大。

kastin

查看过那本书了，那本书给出的结果是对 `n\to \infty`的最大根的，楼主应该在1楼中说明这一条件，不然容易引起误解。

kastin

由于观察到拐点的连线（2楼红点连线）变化规律是s型，只要找出拐点规律，那么在远离对称点的位置处用余切近似求根的时候，可以将余切函数表达式加上一个修正项，即拐点的纵坐标（事实上这对于靠近对称点时求根也可这么做，精度会提高），然后反求出来的根就是精确度非常高的值了。
联想到正切函数或者Logistic函数的反函数可以拟合这种变化，结果如下

1. n = 20;
   
2. datas = {x, f[x, n]} /. NSolve[D[f[x, n], {x, 2}] == 0, x, Reals];
   
3. {datas // ListPlot[#, PlotStyle -> Red] &,
   
4.   Plot[ {Tan[
   
5.       Pi /n (k - n/2)], (-1/a Log[ n/ k - 1] + b /.
   
6.        FindFit[datas, -1/a Log[ n/ k - 1] + b, {a, b}, k])} //
   
7.     Evaluate, {k, 0, n}, PlotStyle -> {Automatic, Orange},
   
8.    PlotLegends -> "Expressions"]} // Show

复制代码

拐点分布

（横轴为`k`，表示的是从左到右第 `k` 个拐点）
上图所示的是问题2的拐点分布，问题1可以同样方法来处理。
修正项为$$-\frac{1}{r}\ln\left(\frac{n}{k}-1\right)+b$$不过这里的待拟合的变量 `r` 和 `b` 与 `k` 的关系是满足类似Logistic曲线规律的（在`n`固定的情况下），如果能找到用 `k` 和 `n` 近似表达`r`和`b`就好了。

不管怎样，上述方法只是对于有拐点的那个分支有作用，对于两端拖尾的那部分曲线，则需要单独处理。

补充内容 (2016-7-7 16:34):
上面的修正项中的 `k` 应该替换为 `k-0.5`（拐点的横坐标非常接近在整数中间）。

282842712474

282842712474 发表于 2016-7-5 10:21
对于1，设$x > n$，则可以用欧拉－麦克劳林求和公式(https://en.wikipedia.org/wiki/Euler%E2%80%93Maclaur ...

不像@kastin 的创造性这么好，我只会按部就班地死算。还是看看这里的方法最远可以走多远。
$$\begin{aligned}a=&-\frac{1}{2x}+\sum_{k=0}^n \frac{1}{x-k}\\
=&-\frac{1}{2x}+\int_0^n \frac{1}{x-k}dk+\frac{1}{2}\left(\frac{1}{x}+\frac{1}{x-n}\right)+\dots\\
=&-\frac{1}{2x}+\ln x-\ln(x-n)+\frac{1}{2}\left(\frac{1}{x}+\frac{1}{x-n}\right)-\frac{1}{12}\left(\frac{1}{x^2}-\frac{1}{(x-n)^2}\right)+\dots
\end{aligned}$$
根据参考答案的形式，我们可以设$A=1-e^{-a}$，即$a=-\ln (1-A)$，得到
$$-\ln(1-A)=-\frac{1}{2x}+\ln x-\ln(x-n)+\frac{1}{2}\left(\frac{1}{x}+\frac{1}{x-n}\right)+\frac{1}{12}\left(\frac{1}{(x-n)^2}-\frac{1}{x^2}\right)+\dots$$
设$x=y/A$，代入上式，在$A=0$处展开得
$$0=\left(\frac{n}{y}+\frac{1}{2 y}-1\right)A +\frac{\left(n^2+n-y^2\right)A^2 }{2 y^2}+\frac{\left(2 n^3+3 n^2+n-2 y^3\right)A^3 }{6 y^3}+\dots$$
只取第一阶，得到
$$0=\left(\frac{n}{y}+\frac{1}{2 y}-1\right)A$$
确实解得
$$y=n+\frac{1}{2}$$
所以首项近似为
$$x=\frac{n+\frac{1}{2}}{1-e^{-a}}$$
逐次逼近，再次设
$$x=\frac{n+\frac{1}{2}}{A}+y$$
代入展开得到
$$\frac{A^2 (-8 n y-4 y-1)}{2 (2 n+1)^2}+\frac{A^3 \left(-24 n^2 y+24 n y^2-24 n y-2 n+12 y^2-1\right)}{3 (2 n+1)^3}+\dots$$
取首项解得
$$y=-\frac{1}{4 (2 n+1)}$$
所以有近似
$$x=\frac{n+\frac{1}{2}}{1-e^{-a}}-\frac{1}{4 (2 n+1)}$$
完整的级数应该是如下的形式
$$x=\frac{n+\frac{1}{2}}{A}-\frac{1}{4 (2 n+1)}+\sum_{k=1}^{\infty} g_k(n) A^k$$
由于$A < 1$，这应该是一个不错的级数。

数学星空

利用下面公式计算：

\(H(k)=\sum_{k=1}^{n} \frac{1}{k}=\ln(n)+\gamma+\frac{1}{2n}-\frac{1}{12n^2}+\frac{1}{120n^4}-\frac{1}{252n^6}+\frac{1}{240n^8}-\frac{1}{132n^{10}}+\dots\)

则有：

\(\frac{1}{2x}+\frac{1}{x-1}+\frac{1}{x-2}+\dots+\frac{1}{x-n}=-\frac{1}{2x}+H(x)-H(x-n-1)\)

\(=-\frac{1}{2x}+\ln(\frac{x}{x-n-1})+(\frac{1}{2x}-\frac{1}{12x^2}+\frac{1}{120x^4}-\frac{1}{252x^6}+\frac{1}{240x^8}-\frac{1}{132x^{10}}+\dots)\)

\(-(\frac{1}{2(n-x-1)}-\frac{1}{12(n-x-1)^2}+\frac{1}{120(n-x-1)^4}-\frac{1}{252(n-x-1)^6}+\frac{1}{240(n-x-1)^8}-\frac{1}{132(n-x-1)^{10}}+\dots)\)

利用上面的方法算出:

记\(A=1-e^{-a}\)

\(x_n=\frac{n+\frac{1}{2}}{A}-\frac{1}{4(2n+1)}+\frac{A(8n^2+8n-1)}{24(2n+1)^2}-\frac{A^2(64n^4+128n^3+60n^2-4n+1)}{48(2n+1)^5}-\frac{A^3(7680n^6+2304n^5+31264n^4+24128n^3+8976n^2+752n+19)}{1440(2n+1)^7}+\dots\)

对于第2问：记\(B=\frac{1-e^{-a}}{1+e^{-a}}\)

\(\frac{1}{x}+\frac{2x}{x^2-1}+\frac{2x}{x^2-2^2}+\dots+\frac{2x}{x^2-n^2}=-\frac{1}{x}+H(x)-H(x-n-1)+H(x+n)-H(x-1)\)

\(x_n=\frac{n+\frac{1}{2}}{B}-\frac{B}{6(2n+1)}-\frac{2B^3(15n^2+15n+1)}{45(2n+1)^3}-\frac{2B^5(1260n^4+2520n^3+861n^2-399n+11)}{945(2n+1)^5}+\dots\)

282842712474

没发现mathematica有自动完成欧拉－麦克劳林求和的函数。

但发现一条途径可以简化过程，充分利用已有结果。我们构造函数（还是对于1）：
$$f(x)=\left(ax+\frac{1}{2}\ln x\right)-\ln \left[x(x-1)\dots (x-n)\right]=\left(ax+\frac{1}{2}\ln x\right)-\ln \frac{\Gamma(x+1)}{\Gamma(x-n)}$$
可见原方程相当于$f'(x)=0$。

因为$\Gamma$函数有斯特灵级数（http://spaces.ac.cn/archives/3731/）可以直接用，因此可以直接用mathematica展开($x\to\infty$)得到
$$\ln \frac{\Gamma(x+1)}{\Gamma(x-n)}=(n+1) \ln x-\frac{n (n+1)}{2 x}-\frac{n \left(2 n^2+3 n+1\right)}{12 x^2}-\frac{n^2 (n+1)^2}{12 x^3}+\dots$$
于是$f'(x)=0$给出：
$$a-\frac{n+\frac{1}{2}}{x}-\frac{n (n+1)}{2 x^2}-\frac{n \left(2 n^2+3 n+1\right)}{6 x^3}-\frac{n^2 (n+1)^2}{4 x^4}+ \dots=0$$
求反演级数得到
$$x=\frac{n+\frac{1}{2}}{a}+\frac{n^2+n}{2 n+1}+\frac{2 a \left(n^4+2 n^3+2 n^2+n\right)}{3 (2 n+1)^3}-\frac{2 a^2 \left(n^4+2 n^3+n^2\right)}{(2 n+1)^5}\dots$$
按照前面的经验，以$a$为参数并不大好，于是可以以$A=1-e^{-a}$为参数，得到
$$x=\frac{n+\frac{1}{2}}{A}-\frac{1}{4 (2 n+1)}+\frac{A \left(8 n^2+8 n-1\right)}{24 (2 n+1)^3}+\frac{A^2 \left(-64 n^4-128 n^3-60 n^2+4 n-1\right)}{48 (2 n+1)^5}+\dots$$

恰好这些东西在mathematica里边都有现成的函数：

1. j = 10;
   
2. G[x_] = (x - 1)!;
   
3. f[x_] = a*x + 1/2*Log[x] -
   
4.    Normal[Simplify[Series[Log[G[x + 1]/G[x - n]], {x, Infinity, j}],
   
5.      Assumptions -> n > 0]];
   
6. g[a_] = Normal[InverseSeries[Series[f'[x], {x, Infinity, j}]]] /.
   
7.    x -> 0;
   
8. Series[g[-Log[1 - A]], {A, 0, j}]

复制代码

282842712474

282842712474 发表于 2016-7-7 07:56
没发现mathematica有自动完成欧拉－麦克劳林求和的函数。

但发现一条途径可以简化过程，充分利用已有结 ...

按照这个思路，可以一鼓作气解决掉2. 感觉2好像更简单，构造$f(x)$
$$f(x)=ax-\ln \prod_{k=-n}^{n} (x+k)$$
那$f'(x)=0$就是原方程了。也就是
$$f(x)=ax-\ln\frac{\Gamma(x+n+1)}{\Gamma(x-n)}$$
修改一下原来的程序

1. j = 10;
   
2. G[x_] = (x - 1)!;
   
3. f[x_] = a*x -
   
4.    Normal[Simplify[
   
5.      Series[Log[G[x + n + 1]/G[x - n]], {x, Infinity, j}],
   
6.      Assumptions -> n > 0]];
   
7. g[a_] = Normal[InverseSeries[Series[f'[x], {x, Infinity, j}]]] /.
   
8.    x -> 0;
   
9. Series[g[-Log[1 - A]], {A, 0, j}]

复制代码

直接给出
$$x=\frac{2 n+1}{A}+\left(-n-\frac{1}{2}\right)-\frac{A}{12 (2 n+1)}-\frac{A^2}{24 (2 n+1)}+\frac{A^3 \left(-120 n^2-120 n-19\right)}{720 (2 n+1)^3}+\dots$$