The Implicit Solution of a Nonlinear Difference Equation
Where the problem comes from
On a math research forum that I follow regularly, there was once a thread discussing the asymptotic solution of the following nonlinear difference equation:
$$a_{n+1}=a_n+\frac{1}{a_n^2},\, a_1=1$$
The original thread is here], and I learned a great deal from it, picking up quite a few new tricks along the way. The main idea was to cube both sides, then set $x_n=a_n^3$, turning the problem into an equivalent recurrence:
$$x_{n+1}=x_n+3+\frac{3}{x_n}+\frac{1}{x_n^2},\,x_1=1$$
From which, through some clever manipulation, an asymptotic expansion could be obtained:
$$x_n = 3n+\ln n+a+\frac{\frac{1}{3}(\ln n+a)-\frac{5}{18}}{n}+\dots$$
I won't repeat the details here; interested readers can go and study the thread above.
However, elegant as this form of the solution is, there are a few points that I find a bit unsatisfying:
1. The solution is an asymptotic series, which means its radius of convergence is actually 0;
2. It's a solution of the form $n^{-k}$, which is hard to evaluate for small $n$, making high-precision computation rather difficult;
3. Of course, the original goal was asymptotic computation, but it doesn't seem necessary to expand so many terms just for asymptotics;
4. It involves a limit value $a$ that is already fairly hard to compute in its own right;
5. The derivation process feels a bit lacking in intuition.
Of course, some of these complaints are nitpicking. But it was precisely these shortcomings that pushed me to look for a solution in a better form, which eventually led to this post. more
The implicit solution
The recurrence
$$x_{n+1}=x_n+3+\frac{3}{x_n}+\frac{1}{x_n^2},\,x_1=1$$
has, as its first-order asymptotic solution, the result obtained by dropping the trailing $\frac{3}{x_n}+\frac{1}{x_n^2}$, giving $x_{n+1}=x_n+3$, and hence $x_n=3n-2$.
There are many possible routes from here — iteration, undetermined coefficients, and so on — all of which can yield asymptotic solutions. But I tried a different route: seeking an implicit solution. First consider:
$$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)$$
Expand the last term to first order around $x_n$:
$$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)$$
Working to order $1/x_n$ and dropping the second-order term, we get
$$x_{n+1}+f(x_{n+1})=x_n+f(x_n)+3+3\left[\frac{1}{x_n}+f'(x_n)\right]$$
If we let $f(x_n)=-\ln x_n$, then the bracketed term vanishes, and we obtain
$$x_{n+1}-\ln x_{n+1}=x_n-\ln x_n+3$$
This equation has precision $\mathcal{O}(x_n^{-2})$; to balance simplicity and precision, it is better to start from $x_2=8$, giving
$$x_n - \ln x_n = 3(n-2)+8-\ln 8$$
This is an implicit (approximate) solution with precision $\mathcal{O}(x_n^{-1})$. It maintains essentially the same precision for almost all $n$. The key point here is that by introducing a new term, the recurrence becomes linear (an arithmetic sequence). This process can be continued: building on the above, we introduce a new function and expand to second order, obtaining:
$$x_{n+1}-\ln x_{n+1}+\frac{5}{6x_{n+1}}=x_n-\ln x_n+\frac{5}{6x_n}+3$$
This has precision $\mathcal{O}(x_n^{-3})$, from which an implicit solution of precision $\mathcal{O}(x_n^{-2})$ can be derived.
Why look for an implicit solution?
Why bother seeking an implicit solution at all? There are, roughly, the following benefits.
Ordinarily, traditional asymptotic solutions are obtained via Taylor series expansion, which itself comes with restrictions (requiring differentiability and a bounded derivative), so it's easy to end up with an asymptotic solution — and even when it isn't strictly asymptotic, the radius of convergence may still be small. An implicit function solution, on the other hand, tends to be more stable and to converge better; even if it is still asymptotic in nature, the rate of divergence is typically much reduced. For instance, $x=\sqrt{1+t}$, when expanded as a power series, has a radius of convergence of only $1$, but can be written in implicit form as $x^2+1=t$, which preserves both precision and simplicity.
In short, anything an explicit expansion can achieve, an implicit function can generally achieve as well; and things an explicit expansion cannot achieve, an implicit function might still manage. So it stands to reason that implicit functions should generally perform better.
Moreover, for the recurrence in this post, the implicit solution is more concise, and does not scatter hard-to-compute instances of $a$ everywhere. Readers might feel that solving a nonlinear equation every time we want a value of $x_n$ sounds complicated. Indeed, if one wishes, one can go straight from the implicit solution to an explicit one by inverting the function — but that just brings us back to an asymptotic series, so there's not much point in doing so.
A recurrence format suited to programming
The two terms computed above were merely an introductory demonstration; since the orders involved were low, the computation wasn't difficult. But to push the computation further — and especially to make it programmable — we need a recurrence format that's easy to reason about, much like the recursive computations used in perturbation methods.
Suppose the $f(x_n)$ we introduce is exact; then clearly, for the exact $f(x)$, we need
$$\frac{3}{x}+\frac{1}{x^2}+f\left(x+3+\frac{3}{x}+\frac{1}{x^2}\right)-f(x)=0$$
At which point the original recurrence becomes
$$x_{n+1}+f(x_{n+1})=x_{n}+f(x_{n})+3$$
which is much easier to solve.
After some analysis, we can artificially introduce a parameter $q$, treat $f(x)$ as a bivariate function $f(x,q)$ of $x,q$, solve for the series solution in $f(x,q)$, and then let $q=1$ to obtain $f(x)=f(x,1)$. The introduced format is as follows:
$$\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$$
The idea behind this construction is: take $x^{-1}$ as the order of the infinitesimal, and let $f^{(n)}(x)$ also be of order $x^{-n}$, so that inside $f()$ we have $x^{-n}\to q^{n+1}x^{-n}$, and outside $f()$ we have $x^{-n} \to q^n x^{-n}$. At this point, software like Mathematica can expand this quickly:
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}]
The result is
$$\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}$$
Setting each order's coefficient of $q$ to zero, we solve successively to get:
$$\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}$$
Hence
$$f(x)=f(x,1)=-\ln x+\frac{5}{6x}+\frac{2}{3 x^2}+\dots$$
That is,
$$\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}$$
giving
$$\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}$$
which has precision $\mathcal{O}(x_n^{-3})$.
With a slight change to this approach, we can write a Mathematica program to automate the computation:
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
By changing the 5 in {n, 0, 5}, one can raise or lower the precision. The code above gives:
$$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$$
Asymptotic results and the limit value
This gives us
$$x_n+f(x_n)=3(n-k)+x_k+f(x_k),\,\text{with}x_k\text{with initial value}$$
Let
$$a=x_k+f(x_k)-3k+\ln 3$$
which gives
$$x_n + f(x_n)=3n+a-\ln 3$$
To obtain an explicit asymptotic expression, we can build an iteration on top of the above:
$$x^{(k+1)}_n=3n+a-\ln 3 -f(x^{(k)}_n)$$
Here $x^{(k)}_n$ denotes the $k$-th order approximation of $x_n$. Starting from $x^{(0)}_n=3n+a-\ln 3$ as the initial value and iterating, then expanding, we get
$$\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}$$
This matches the asymptotic expression given by mathe on the research forum.
The Mathematica code for the iteration is (continuing from before):
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}]
From this we discover that $a$ here turns out to be exactly the limit value $a$:
$$a=\lim_{n\to\infty} (x_n - 3n - \ln n)$$
This simultaneously gives us the asymptotic expression for $a$:
$$x_k+f(x_k)-3k+\ln 3$$
Using the value of $x_{10^9}$ computed by wayne, and substituting it into the expression above, we obtain
$$a=1.1352558473155037141953943477479\dots$$
Translated automatically with claude-sonnet-5; all equations are reproduced verbatim from the source. Copyright remains with the original author.