Finding solution of $\frac{e^x}{2} = 1 + x + \frac{x^2}{2!} + \dots + \frac{x^k}{k!}$

Question

I am looking to find a numerical solution to the equation:

$$\frac{e^x}{2} = 1 + x + \frac{x^2}{2!} + \dots + \frac{x^{k-1}}{(k-1)!}$$

where $k$ is of the order of $10^7$, and the solution is accurate to at least 2-3 decimals. It is not hard to show that the equation has exactly one positive solution (just differentiate and use induction).

Roughly, it measures how good of an approximation is the Taylor expansion of $e^x$ upto $k$ terms. Thus, as $k$ increases, we can expect $x$ to increase as well, as the accuracy of the expansion increases farther from the origin. We can get an lower bound on $x$ by the Taylor's theorem on $e^x$, as the remainder here is $\frac{e^x}{2}$, so $\frac{e^x}{2} \leq \frac{e^x}{k!}x^k$ which implies $\sqrt[k]{\frac{k!}{2}} \leq x$. Thus by Stirling's formula, we see that the lower bound is almost $\frac{k}{e}$.

My first idea was to use binary search, by simply computing the entire function and checking its sign. I coded it up in Python, but unfortunately, it only works for values of $k$ upto about $1000$. Beyond that, the values get too big to fit in its numeric data type. I tried many other ways, like Newton's method which did not work out, amongst other things. The main problem here is I can't find a way to avoid computing the entire function (which will overflow). I have tried hard to solve it, but couldn't do it so hopefully you all can help.

Answer 1

This is the exponential sum function. Letting $k=n+1$ and using identities with the incomplete Gamma function and the regularized Gamma function, the problem can be rewritten as solving:

$$Q(n+1,x)-\frac12=0$$

The problem is then tractable, as the main problem is numerically evaluating the function due to a significant amount of cancellation and the fact that $e^x$ simply becomes $\infty$ near the root in double precision as you have said.

We can now evaluate this using WolframAlpha in the above form. In only a few secant iterations, we find that the desired root for $n=10^7$ is then

$$x=10000000.666666668\dots$$

---

It is interesting to note that the root seems to remain between $n$ and $n+1$, which makes sense when working out the Stirling approximation with the Taylor remainder.

Answer 2

From a numerical point of view, I should write the problem as : find the zero of $$f(\epsilon)=\log \left(\frac{\Gamma (k+1,k+\frac 23+\epsilon}{k!}\right)+\log (2)$$ A quite detailed numerical analysis shows that $$\epsilon=\frac{8}{405\, k}-\frac{15}{1196 \,k^2}+\frac{7}{1170\, k^3}-\frac{1}{722 \, k^4}$$ is a quite good approximation (all parameters being highly significant). So, as an approximation $$\color{blue}{x=k+\frac 23+\frac{8}{405\, k}-\frac{15}{1196 \,k^2}+\frac{7}{1170\, k^3}-\frac{1}{722 \, k^4}}$$

Applied to the case $k=10^7$ used by @Simply Beautiful Art, the above formula gives $$x=\color{red}{1.000000066666666864197518322}39\times 10^7$$ while the exact solution is $$x=1.00000006666666686419751832256\times 10^7$$

Edit

From a formal point of view, the first iterate of Newton method gives $$x_1=x_0+\frac{e^{x_0}}{ x_0^{k}}\, \Gamma (k+1,x_0)\,\log \left(\frac{\Gamma (k+1,x_0)}{k!}\right)\qquad \text{with}\qquad x_0=k+\frac 23$$

Computing $\epsilon$ for $5 \leq k \leq 10000$ (step size = $5$), a polynomial regression gives (after rationalization of the coefficients) $$\epsilon=\frac{8}{405\, k}-\frac{15}{1198\, k^2}+\frac{8}{1299 \,k^3}-\frac{2}{1007\, k^4}$$ $$\begin{array}{clclclclc} \text{} & \text{Estimate} & \text{Standard Error} & \text{Confidence Interval} \\ a & +0.019753 & 2.76 \times 10^{-8} & \{+0.019753,+0.019753\} \\ b & -0.012521 & 1.42 \times 10^{-6} & \{-0.012524,-0.012518\} \\ c & +0.006159 & 1.81 \times 10^{-5} & \{+0.006124,+0.006194\} \\ d & -0.001987 & 5.87 \times 10^{-5} & \{-0.002102,-0.001871\} \\ \end{array}$$

We could even do better using Halley method which would give $$x_2=x_0+\frac{2 f(x_0) f'(x_0)}{f(x_0) f''(x_0)-2 f'(x_0)^2}$$ $$f'(x_0)=-\frac{e^{-x_0} \left(x_0\right)^k}{\Gamma \left(k+1,x_0\right)}\qquad f''(x_0)=\frac{e^{-2 x_0} x_0^{k-1} \left(2 e^{x_0} \Gamma (k+1,x_0)-3 x_0^{k+1}\right) } {3 \Gamma (k+1,x_0)^2 }$$

Update

Thanks to the link @Semiclassical provided, the already existing solution is $$\epsilon=\frac{8}{405\, k}-\frac{64}{5103\, k^2}+\frac{2944}{492075 \, k^3}+O\left(\frac{1}{k^4}\right)$$

Continuing my work, it seems that the next term could be close to $$- \frac{149}{113289\,k^4}$$

Answer 3

This answer says $$ \left[\sum_{k=0}^n\frac{x^k}{k!}-\frac12e^x\,\right]\overset{\substack{x=n\\[2pt]\\}}=\frac23\frac{e^n}{\sqrt{2\pi n}}+O\!\left(\frac{e^n}n\right)\tag1 $$ Note that $$ \begin{align} \frac{\mathrm{d}^m}{\mathrm{d}x^m}\left[\sum_{k=0}^n\frac{x^k}{k!}-\frac12e^x\,\right] &=\left[\sum_{k=0}^n\frac{x^k}{k!}-\frac12e^x\,\right]-\sum_{k=n-m+1}^n\frac{x^k}{k!}\tag{2a}\\[3pt] &\overset{\substack{x=n\\[2pt]\\}}=\frac23\frac{e^n}{\sqrt{2\pi n}}-m\frac{e^n}{\sqrt{2\pi n}}+O\!\left(\frac{e^n}n\right)\tag{2b} \end{align} $$ Applying Taylor's Series, we get $$ \begin{align} \sum_{k=0}^n\frac{(n+x)^k}{k!}-\frac12e^{n+x} &=\frac{e^n}{\sqrt{2\pi n}}\sum_{m=0}^\infty\left(\frac23-m\right)\frac{x^m}{m!}+O\!\left(\frac{e^n}n\right)\tag{3a}\\ &=\left(\frac23-x\right)\frac{e^{n+x}}{\sqrt{2\pi n}}+O\!\left(\frac{e^n}n\right)\tag{3b} \end{align} $$ Thus, as $n\to\infty$, it looks as if $x=n+\frac23$ is pretty close.

---

Improvement

As noted in a comment to the answer cited above, $$ \sum_{k=0}^n\frac{x^k}{k!}-\frac12e^x\overset{\substack{x=n\\[2pt]\\}}= \frac{e^n}{\sqrt{2\pi n}}\left(\frac23-\frac{23}{270n}+\frac{23}{3024n^2}+\frac{259}{77760n^3}+O\!\left(\frac1{n^4}\right)\right)\tag4 $$ $(4)$ extends $(1)$ and part of $\text{(2a)}$. To extend the other part of $\text{(2a)}$, we use Stirling Numbers of the First Kind: $$\newcommand{\stirone}[2]{\left[{#1}\atop{#2}\right]} \begin{align} &\sum_{k=0}^{m-1}\frac{x^{n-k}}{(n-k)!}\tag{5a}\\ &\overset{\substack{x=n\\[2pt]\\}}=\frac{n^n}{n!}\sum_{k=0}^{m-1}\prod_{j=0}^{k-1}\left(1-\frac{j}{n}\right)\tag{5b}\\ &=\frac{n^n}{n!}\sum_{k=0}^{m-1}\sum_{j=0}^k\stirone{k}{k-j}\left(-\frac1n\right)^j\tag{5c}\\ &=\frac{n^n}{n!}\sum_{j=0}^{m-1}\sum_{k=j}^{m-1}\stirone{k}{k-j}\left(-\frac1n\right)^j\tag{5d}\\ &=\frac{n^n}{n!}\left(\binom{m}{1}-\frac{\binom{m}{3}}n+\frac{3\binom{m}{5}+2\binom{m}{4}}{n^2}-\frac{15\binom{m}{7}+20\binom{m}{6}+6\binom{m}{5}}{n^3}+O\!\left(\frac1{n^4}\right)\right)\tag{5e} \end{align} $$ Explanation:
$\text{(5a)}$: the right side of the difference in $\text{(2a)}$
$\text{(5b)}$: expand as a product
$\text{(5c)}$: apply Stirling Numbers of the First Kind
$\text{(5d)}$: switch order of summation
$\text{(5e)}$: use Stirling Number identities to sum the terms for $0\le j\le3$

Summing $\text{(5e)}$ against $\frac{x^m}{m!}$ gives $$ \frac{n^n}{n!}e^x\!\left(x-\frac{x^3}{6n}+\frac{3x^5+10x^4}{120n^2}-\frac{15x^7+140x^6+252x^5}{5040n^3}+O\!\left(\frac1{n^4}\right)\right)\tag6 $$ Taking the reciprocal of Stirling's Formula gives $$ \frac{n^n}{n!}=\frac{e^n}{\sqrt{2\pi n}}\left(1-\frac1{12n}+\frac1{288n^2}+\frac{139}{51840n^3}+O\!\left(\frac1{n^4}\right)\right)\tag7 $$ Combining $(6)$ and $(7)$, subtracting from $(4)$ summed against $\frac{x^m}{m!}$, gives, via Taylor's Theorem, $$ \begin{align} &\sum_{k=0}^n\frac{(n+x)^k}{k!}-\frac12e^{n+x}\\ &=\frac{e^{n+x}}{\sqrt{2\pi n}}\left(\vphantom{\frac{x^2}{n^2}}\right.\frac{2-3x}3+\frac{90x^3+45x-46}{540n}\\ &-\frac{756x^5+2520x^4+420x^3+105x-230}{30240n^2}\\[3pt] &-\frac{3240x^7+30240x^6+52164x^5-7560x^4-630x^3+2919x-3626}{1088640n^3}\left.\vphantom{\frac{x^2}{n^2}}+O\!\left(\frac1{n^4}\right)\right)\tag8 \end{align} $$