問題來源#

筆者經常學習的數學研發論壇曾有一帖討論下述非線性差分方程的漸近求解：

$$a_{n+1}=a_n+\frac{1}{a_n^2},\, a_1=1$$

原帖子在這里，從這帖子中我獲益良多，學習到了很多新技巧。主要思路是通過將兩邊立方，然后設$x_n=a_n^3$，變為等價的遞推問題：

$$x_{n+1}=x_n+3+\frac{3}{x_n}+\frac{1}{x_n^2},\,x_1=1$$

然后可以通過巧妙的技巧得到漸近展開式：

$$x_n = 3n+\ln n+a+\frac{\frac{1}{3}(\ln n+a)-\frac{5}{18}}{n}+\dots$$

具體過程就不提了，讀者可以自行到上述帖子學習。

然而，這種形式的解雖然精妙，但存在一些筆者不是很滿意的地方：1、解是漸近的級數，這就意味著實際上收斂半徑為0；

2、是$n^{-k}$形式的解，對于較小的$n$難以計算，這都使得高精度計算變得比較困難；

3、當然，題目本來的目的是漸近計算，但是漸近分析似乎又沒有必要展開那么多項；

4、里邊帶有了一個本來就比較難計算的極限值$a$；

5、求解過程似乎稍欠直觀。

當然，上面這些缺點，有些是雞蛋里挑骨頭的。不過，也正是這些缺點，促使我尋找更好的形式的解，最終導致了這篇文章。

隱式解#

遞推

$$x_{n+1}=x_n+3+\frac{3}{x_n}+\frac{1}{x_n^2},\,x_1=1$$

的第一級漸近解是通過略去后面的$\frac{3}{x_n}+\frac{1}{x_n^2}$得出的，得到$x_{n+1}=x_n+3$，所以$x_n=3n-2$。

接下來的思路有很多，比如迭代、待定系數等，都可以求得漸近解，但是我嘗試了另外一條途徑：求隱式的解。首先我們考慮：

$$x_{n+1}+f(x_{n+1})=x_n+3+\frac{3}{x_n}+\frac{1}{x_n^2}+f\left(x_n+3+\frac{3}{x_n}+\frac{1}{x_n^2}\right)$$

將最后一項在$x_n$處展開至1階：

$$x_{n+1}+f(x_{n+1})=x_n+3+\frac{3}{x_n}+\frac{1}{x_n^2}+f(x_n)+f'(x_n)\left(3+\frac{3}{x_n}+\frac{1}{x_n^2}\right)$$

以$1/x_n$為階，略去二階項，得到

$$x_{n+1}+f(x_{n+1})=x_n+f(x_n)+3+3\left[\frac{1}{x_n}+f'(x_n)\right]$$

如果讓$f(x_n)=-\ln x_n$，那么方括號就為0，于是得到

$$x_{n+1}-\ln x_{n+1}=x_n-\ln x_n+3$$

這個等式的精度是$\mathscr{O}(x_n^{-2})$，為了兼顧簡潔與精度，從$x_2=8$出發比較好，這樣

$$x_n - \ln x_n = 3(n-2)+8-\ln 8$$

這就是一個隱式(近似)解，精度是$\mathscr{O}(x_n^{-1})$。它幾乎對于所有的$n$都保持同樣的精度。這里的要點就是引入新的項，使得它變成了線性的遞推(等差數列)。這個過程可以繼續下去，在前述的基礎上引入新的函數，繼續將它展開到二階，算得到：

$$x_{n+1}-\ln x_{n+1}+\frac{5}{6x_{n+1}}=x_n-\ln x_n+\frac{5}{6x_n}+3$$

它具有$\mathscr{O}(x_n^{-3})$的精度，由此可以解得精度為$\mathscr{O}(x_n^{-2})$的隱式解。

為什么要隱式解？#

為什么要尋找隱式解？大概有如下的好處。

一般情況下，傳統漸近解是通過泰勒級數展開得到的，泰勒級數展開本身就存在限制(要求可導且導函數有限)，因此容易出現漸近解，即使不漸近，收斂半徑也可能較小。而隱函數的解通常來說比較穩定，收斂性比較好，雖然可能還是漸近的，但是發散速度會降低不少。比如$x=\sqrt{1+t}$，展開為冪級數，收斂半徑僅為$1$，但是可以表示為隱函數形式$x^2+1=t$，保持了精度和簡潔。

簡而言之，如果顯式展開式能做到的，基本上隱函數也能做到；而顯式展開式不能做到的，隱函數也可能做到。因此顯然隱函數性能會更好。

此外，對于本文的遞推來說，隱函數的解更加簡潔，不到處出現難以計算的$a$。可能讀者覺得，每求一次$x_n$都需要求解一次非線性方程，比較復雜。事實上，如果愿意的話，可以直接通過求反函數來從隱式解獲得顯式解，但是這又回到了漸近級數了，就沒有什么必要了。

利于編程的遞推格式#

前面所計算的兩項，僅僅是介紹性的演示，由于階數不高，計算也不困難，但是為了進一步計算下去，尤其是為了編程計算，需要構造便于理解的遞推格式，就好比攝動法的遞推計算。

假設引入的$f(x_n)$是精確的，那么顯然，對于精確的$f(x)$，需要滿足

$$\frac{3}{x}+\frac{1}{x^2}+f\left(x+3+\frac{3}{x}+\frac{1}{x^2}\right)-f(x)=0$$

這時候原來的遞推就變成了

$$x_{n+1}+f(x_{n+1})=x_{n}+f(x_{n})+3$$

求解就容易多了。

經過分析，可以通過人工引入參數$q$，并且將$f(x)$當作$x,q$的二元函數$f(x,q)$，求得$f(x,q)$的級數解，然后讓$q=1$，得到$f(x)=f(x,1)$。引入的格式如下：

$$\frac{3q}{x}+\frac{q^2}{x^2}+f\left(x+3q+\frac{3q^2}{x}+\frac{q^3}{x^2},q\right)-f(x,q)=0$$

這樣引入的思路是：以$x^{-1}$為無窮小階，且讓$f^{(n)}(x)$也具有$x^{-n}$的階，所以，在$f()$內$x^{-n}\to q^{n+1}x^{-n}$，在$f()$外有$x^{-n} \to q^n x^{-n}$。這時候，可以用Mathematica之類的軟件，很快展開它：Series[f[x + 3\cdot q + 3\cdot q^2/x + q^3/x^2, q] - f[x, q] + 3\cdot q/x +

q^2/x^2, {q, 0, 5}]

結果是

$$\begin{aligned}&q \left(3 f^{(1,0)}(x,0)+\frac{3}{x}\right)\\

+&q^2 \left(\frac{3 f^{(1,0)}(x,0)}{x}+3 f^{(1,1)}(x,0)+\frac{9}{2} f^{(2,0)}(x,0)+\frac{1}{x^2}\right)\\

+&q^3 \left(\frac{f^{(1,0)}(x,0)}{x^2}+\frac{3 f^{(1,1)}(x,0)}{x}+\frac{3}{2} f^{(1,2)}(x,0)\right.\\

&\qquad\qquad\left.+\frac{9 f^{(2,0)}(x,0)}{x}+\frac{9}{2} f^{(2,1)}(x,0)+\frac{9}{2} f^{(3,0)}(x,0)\right)\\

+&\dots

\end{aligned}$$

讓$q$的各階系數為0，逐次解得：

$$\begin{aligned}&f^{(0,0)}(x,0)=-\ln x\\

&f^{(0,1)}(x,0)=\frac{5}{6x}\\

&f^{(0,2)}(x,0)=\frac{4}{3 x^2}\\

&\dots

\end{aligned}$$

因此

$$f(x)=f(x,1)=-\ln x+\frac{5}{6x}+\frac{2}{3 x^2}+\dots$$

也就是

$$\begin{aligned}&x_{n+1}-\ln x_{n+1}+\frac{5}{6x_{n+1}}+\frac{2}{3 x_{n+1}^2}\\

=&x_{n}-\ln x_{n}+\frac{5}{6x_{n}}+\frac{2}{3 x_{n}^2}+3\end{aligned}$$

得到

$$\begin{aligned}&x_{n}-\ln x_{n}+\frac{5}{6x_{n}}+\frac{2}{3 x_{n}^2}\\

=&3(n-2)+8-\ln 8+\frac{5}{6\times 8}+\frac{2}{3\times 8^2}\end{aligned}$$

這具有$\mathscr{O}(x_n^{-3})$的精度。

將思路稍加變動，就可以用Mathematica寫出自動計算程序：g[x_] = 0;

ff[x_] = 3\cdot 1/x + 1/x^2

Do[e = q^n;

g[x_] = g[x] +

q^n\cdot Integrate[-D[

f[x + 3\cdot q\cdot e + ff[x/q]\cdot e\cdot q, q] - f[x, q] +

g[x + 3\cdot q + ff[x/q]\cdot q] - g[x] + ff[x/q], {q,

n + 1}] /. {q -> 0, f -> 0} // Simplify, x]/

3/(n + 1)!;, {n, 0, 5}]

h[x_] = g[x] /. {q -> 1} // Expand

通過修改{n, 0, 5}中的5可以增減精度。上述代碼給出：

$$f(x)=-\ln x+\frac{5}{6 x}+\frac{2}{3 x^2}+\frac{77}{108 x^3}+\frac{133}{240 x^4}-\frac{2669}{5400 x^5}＋\dots$$

漸近結果與極限值#

此時解得

$$x_n+f(x_n)=3(n-k)+x_k+f(x_k),\,\text{以}x_k\text{為初值的情況下}$$

設

$$a=x_k+f(x_k)-3k+\ln 3$$

得到

$$x_n + f(x_n)=3n+a-\ln 3$$

如果要得到顯式的漸近表達式，可以以上式為基礎構建迭代式：

$$x^{(k+1)}_n=3n+a-\ln 3 -f(x^{(k)}_n)$$

這里的$x^{(k)}_n$指的是$x_n$的第$k$級近似。以$x^{(0)}_n=3n+a-\ln 3$為初值，進行迭代，并展開，得到

$$\begin{aligned}x_n=&3n+a+\ln n+\frac{6 a+6 \ln n-5}{18 n}+\\

&\frac{-3 a^2-6 a \ln n+11 a-3 \ln ^2 n+11 \ln n-9}{54 n^2}+\dots\end{aligned}$$

這跟研發論壇上mathe給出的漸近式是一樣的。

迭代的Mathematica代碼為(接前)：p[n_] = 3\cdot n + a - Log[3];

Do[p[n_] =

Normal[Series[-h[p[n]] + 3\cdot n + a - Log[3], {n, Infinity, 5}]], {i,

0, 5}]

Series[p[n], {n, Infinity, 5}]

由此發現，這里的$a$正好是極限值$a$：

$$a=\lim_{n\to\infty} (x_n - 3n - \ln n)$$

這也正好同時給出了求$a$的漸近表達式：

$$x_k+f(x_k)-3k+\ln 3$$

利用wayne算出的$x_{10^9}$，代入上式，算得

$$a=1.1352558473155037141953943477479\dots$$

更詳細的轉載事宜請參考：《科學空間FAQ》

如果您還有什么疑惑或建議，歡迎在下方評論區繼續討論。

如果您覺得本文還不錯，歡迎分享/打賞本文。打賞并非要從中獲得收益，而是希望知道科學空間獲得了多少讀者的真心關注。當然，如果你無視它，也不會影響你的閱讀。再次表示歡迎和感謝！

打賞

微信打賞

支付寶打賞

因為網站后臺對打賞并無記錄，因此歡迎在打賞時候備注留言。你還可以點擊這里或在下方評論區留言來告知你的建議或需求。

如果您需要引用本文，請參考：

蘇劍林. (Apr. 09, 2016). 《一個非線性差分方程的隱函數解 》[Blog post]. Retrieved from https://spaces.ac.cn/archives/3696

總結

以上是生活随笔為你收集整理的php类中引函数变量,一个非线性差分方程的隐函数解的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。