I computed the integral of $\sqrt{1+e^x}$ by hand and got
$$2\sqrt{1+e^x} + \ln\left(\sqrt{1+e^x} – 1\right) – \ln\left(\sqrt{1+e^x} + 1\right) + C,$$
or
$$2\sqrt{1+e^x} + \ln\left(\frac{\left(\sqrt{1+e^x} – 1\right)^2}{e^x}\right) + C.$$
But when I plot the graph on Desmos, this strange thing happened:
Can anyone explain this discontinuity phenomenon I have seen? And, is there a better form of the function that the one I have?
Best Answer
The problem is very likely due to the finite precision of the arithmetic that's used (presumably IEEE-754 double), together with bad conditioning due to cancellation:
In $\ln(\sqrt{1+e^x}-1)$, when $x$ is small ($x\to-\infty$) then $e^x\to 0$ and the square-root evaluates to $\approx 1+e^x/2$. The magnitude of that term is $\approx 1$, and when you subtract $1$ then what you get is $e^x/2$ but only with the precision in which $1$ was represented. This is because the precision of $1+e^x/2$ is basically the precision of $1$ when $e^x\to0$.
When you coalesced the two $\ln$ terms, you got
$$\begin{align} \ell(x) = \ln\left(\sqrt{1+e^x} - 1\right) - \ln\left(\sqrt{1+e^x} + 1\right) &= \ln\left(\frac{\sqrt{1+e^x} - 1}{\sqrt{1+e^x} + 1}\right) \\ &= \ln\left(\frac{(\sqrt{1+e^x} - 1)^2}{e^x}\right) \end{align}$$ which made the problem of cancellation in $\sqrt{1+e^x}-1$ even worse because you squared that term to get a "nice" denominator. Provided $e^x$ is small, squaring that term means doubling $x$, thus $x$ gets even more negative which amplifies the problem.
Going for a nice numerator, however, gives:
$$\begin{align} \ell(x) &= \ln\left(\frac{e^x}{(\sqrt{1+e^x} + 1)^2}\right) \\ &= x - 2\ln(\sqrt{1+e^x}+1) \\\tag 2 \end{align}$$ The representation in $(2)$ mitigated the problem of cancellation, and it's still a nice one. The nasty $\sqrt{1+e^x}-1$ has disappeared altogether.
Note: I cross-checked my result against Desmos. The artifact is around $x=-35$:
https://www.desmos.com/calculator/wxj6kiofa4
https://i.stack.imgur.com/R45tz.png
Where to expect the breakdown
The smallest value greater than $1$ that can be represented by IEEE-754 double is $1+\varepsilon$ with $\varepsilon=2^{-52}$. Determining $x$ such that $\sqrt{1+e^x}-1 = \varepsilon$ means $$1+e^x = (1+\varepsilon)^2 \approx 1+2\varepsilon$$ $$x \approx \ln (2\varepsilon) \approx -51\ln2\approx-35.4$$
And you actually observed a breakdown at $x\approx-35$.