数学研发网设为首页收藏本站

数学研发论坛

 找回密码
 欢迎注册
查看: 645|回复: 27

[原创] 一个代数方程根的级数解问题

[复制链接]
发表于 2016-7-4 20:47:55 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?欢迎注册

x
求下列关于\(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}}\)

我们现在需要做的是给出更精确的级数解。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-7-5 10:21:14 | 显示全部楼层
本帖最后由 282842712474 于 2016-7-5 10:33 编辑

对于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$很小的时候。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-7-5 10:26:37 | 显示全部楼层
本帖最后由 282842712474 于 2016-7-5 10:29 编辑

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

点评

@数学星空,n+1/2让我想到了Gamma函数,有几个公式中自变量恰好就是n+1/2. 见https://en.wikipedia.org/wiki/Gamma_function  发表于 2016-7-5 18:38
@数学星空 这不难理解,多考虑后面一项即可,可是我不是很喜欢先求和后求反函数的迂回的方法,所以在想新的思路  发表于 2016-7-5 11:23
问题的答案就在于如何理解分子\(\frac{1}{2}+n\),此问题来源于《数学分析中的定理和问题》  发表于 2016-7-5 10:48
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-7-5 20:24:37 | 显示全部楼层
本帖最后由 kastin 于 2016-7-5 20:38 编辑

这里有个逆级数链接 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

问题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

问题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` 的取值影响更大。

评分

参与人数 1贡献 +6 经验 +6 鲜花 +6 收起 理由
数学星空 + 6 + 6 + 6 赞一个!

查看全部评分

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-7-6 14:38:11 | 显示全部楼层
本帖最后由 kastin 于 2016-7-6 14:41 编辑

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

点评

原题只针对n→∞,本题需要考虑各个根存在的区间考虑近似根公式哈,难度很大!  发表于 2016-7-6 15:27
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-7-6 18:37:28 | 显示全部楼层
由于观察到拐点的连线(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`(拐点的横坐标非常接近在整数中间)。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-7-6 20:53:27 | 显示全部楼层
本帖最后由 282842712474 于 2016-7-6 20:56 编辑
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$,这应该是一个不错的级数。

点评

你的形式解中\(g_k(n)\)求解也不太容易  发表于 2016-7-7 07:19

评分

参与人数 1金币 +8 贡献 +6 经验 +6 鲜花 +6 收起 理由
数学星空 + 8 + 6 + 6 + 6 方法很有启发性,可以多算几项,然后数值计.

查看全部评分

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2016-7-6 22:33:36 | 显示全部楼层
利用下面公式计算:

\(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\)

点评

不错,用$B$来做参数~我忽略了这一点了。  发表于 2016-7-7 19:50
重新计算了一下,的确也可以算出答案  发表于 2016-7-7 19:16
但是仔细想了一下,这个思路估计可以,你展开式写错了。“则有”后面式子的等号右边,第二个圆括号里边,自变量应该是x不是n  发表于 2016-7-7 08:03
嗯,你说的有道理  发表于 2016-7-7 07:17
不能这样用吧?调和级数要求n是整数,这里$x$一般是无理数  发表于 2016-7-7 06:53
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-7-7 07:56:04 | 显示全部楼层
本帖最后由 282842712474 于 2016-7-7 08:20 编辑

没发现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}]
复制代码

评分

参与人数 1金币 +12 贡献 +12 经验 +12 鲜花 +12 收起 理由
数学星空 + 12 + 12 + 12 + 12 很给力!

查看全部评分

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-7-7 08:13:52 | 显示全部楼层
本帖最后由 282842712474 于 2016-7-7 08:24 编辑
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$$
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|Archiver|数学研发网 ( 苏ICP备07505100号  

GMT+8, 2017-8-20 08:15 , Processed in 0.274855 second(s), 31 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.

快速回复 返回顶部 返回列表