What does one do when Wolframalpha fails to solve an integral? I'm looking for an analytic solution to this integral:
$$
\int\sqrt{\left(\frac{2 d}{(1 + 4 x ^ 2)^{\frac{3}{2}}} + 1\right) \frac{2}{\sqrt{1 + 4 x ^ 2}}}dx, d > -0.5
$$
Wolframalpha fails to find a solution to this, and using quadratures for evaluation is too slow for my use case since I need this for a real-time application.
What does one do in such cases? My intuition is to try to find a function that approximates the integral decently enough on a graphic calculator, but it's quite challenging to start from scratch.
The challenging part of this integral is the $+ 1$ in the first term. Without it, the solution is trivial.
Update:
The closest approximation I found is this:
$$
\frac{6 x \sqrt{\frac{2 d}{\sqrt{1 + 4 x ^ 2}} + 1}}{\sqrt[4]{9 ^ 2 x ^ 2 + 8 ^ 2} + 2}
$$
How could I improve the accuracy of this?
 A: First of all:
$$\sqrt{\frac{2d}{(4x^2+1)^\frac32}+1}=\sum_{n=0}^\infty\binom{\frac12}n\left(\frac{2d}{(4x^2+1)^\frac32}\right)^n= \sum_{n=0}^\infty\binom{\frac12}n(2d)^n(4x^2+1)^{-\frac{3n}2}$$
Therefore:
$$\int\sqrt{\left(\frac{2 d}{(1 + 4 x ^ 2)^{\frac{3}{2}}} + 1\right)} \frac{\sqrt2}{\sqrt[4]{1 + 4 x ^ 2}} dx= \sqrt2\sum_{n=0}^\infty\binom{\frac12}n(2d)^n\int\frac1{\sqrt[4]{4 x ^ 2+1}}(4x^2+1)^{-\frac{3n}2}dx$$
we integrate with the incomplete beta function $\text B_z(a,b)$:
$$\int\sqrt{\frac2{\sqrt{4x^2+1}}\left(\frac{2 d}{(1 + 4 x ^ 2)^{\frac{3}{2}}} + 1\right)}dx=-\frac{i\text{sgn}(x)}{2^\frac32}\sum_{n=0}^\infty\binom{\frac12}n (2d)^n \text B_{-4x^2}\left(\frac12,\frac34-\frac{3n}2\right)$$
which works when comparing both test links’ results, where the $5$th term in the example is on the order of $10^{-8}$ with potentially faster decaying terms as $n\to\infty$. It also works for the binomial series radius of convergence $2|d|<(4x^2+1)^\frac32\iff d\ge\frac12,x>\pm\frac12\sqrt{(2d)^\frac23-1}$:

Now you have an analytic solution as requested
A: This is what we have done in a context where we had to compute the solution of an equation system in $\mathbb{R}^4$, with time constraints like you have.

*

*Numerically (and offline), compute the integral for $N$ values of $d$, and the integration bounds if they vary.

*Find an approximation of the function which minimizes the maximum error on this set of $N$ values. You can use polynomials, linear combinations of positive and negative powers of variables, rational functions,... whatever seems more fitting and will not be too long to compute in your real-time application. You may probably find some python package to help finding the set of parameters that minimize the maximum error. There is a theory of Padé approximant that may be relevant too, although I have not used it.

*Check the validity of the approximation on another set of points.

