How to create a Pade approximation for a difficult function with a divergent Taylor series?

Question

I've been trying to create a good approximation for this function:

$$h\left(a\right)\equiv-\int_{-\infty}^{\infty}p\left(x,a\right)\ln\left(p\left(x,a\right)\right)dx$$ where
$$p\left(x,a\right)\equiv\frac{1}{2\sqrt{2\pi}}\left(e^{-(x+a)^{2}/2}+e^{-(x-a)^{2}/2}\right)$$

This is just the entropy of a bimodal distribution, expanding in the distance of the modes from the center. I haven't been able to figure out the domain of convergence of the Taylor series, but it appears to be quite small:

$$h\left(a\right)=\frac{1+\ln\left(2\pi\right)}{2}+\frac{1}{2}a^{2}-\frac{1}{4}a^{4}+\frac{1}{6}a^{6}-\frac{5}{24}a^{8}+\frac{13}{30}a^{10}-\frac{227}{180}a^{12}+\frac{2957}{630}a^{14}\\-\frac{21425}{1008}a^{16}+\frac{642853}{5670}a^{18}+\ldots$$

I thought maybe a Pade approximation might be better than Taylor, but I always see those calculated from the Taylor Series, which won't work well here. So what can I do?

One thing I tried began with noting that $h(a)$ has a rather Gaussian shape, and so $\ln\left(\frac{h(\infty)-h(a)}{h(\infty)-h(0)}\right)$ is nearly a parabola. But the $a^2$ term of that expansion wasn't even a good fit to the parabola. But maybe a Padme approximation would fit this function better than it would to $h(a)$?

I'm really at a loss. There's something about this function that makes it difficult to even approximate well.

Answer 1

For small values of $a$, the Pade approximants $P^N_M$ calculated from the terms of your power series give a good approximation. Here is a plot of the (log of) ratio of some diagonal Pade approximants/numeric:

This is more of a back-of-envelope calculation: you should play around with it and also look at the $P^N_{N+1}$.

For large values of $a$, split the integral up

$$\tag{1} 2 \sqrt{2\pi} \ h(a)=\int dx \ e^{-(x+a)^2/2}\ln(p) \ +\int dx \ e^{-(x-a)^2/2}\ln(p) $$

Looking at just the first term on the RHS for now, change integration variables $u=x+a$, and expand $\ln (p)$ near the maxima at $u=0$ to get

$$\tag{2} \int du \ e^{-u^2/2}\left[p_0+u^2p_2+\cdots \right] $$

Only terms even in $u$ will survive the integral. The $p_i(a)$ are coefficients of expanding $\ln(p)$. Notice that in the second term (1) we may change variables $u=x-a$ to get an expression identical to (2), so we have

$$\tag{3} h(a)\sim \frac{1}{\sqrt{2\pi}}\int du \ e^{-u^2/2}\left[p_0+u^2p_2+\cdots \right] \qquad, \qquad a \to \infty $$

Term by term integration (thanks Mathematica!) yields

$$\tag{4} h(a)\sim \frac{\ln(8\pi)}{2}-\ln(1+e^{-2a^2})+\frac{1}{2}(1-a^2\operatorname{sech}^2(a^2)) \qquad ,\qquad a\to \infty $$

Here is a plot of the ratio approx/numeric. Note the vertical scale:

It turns out the approximation (4) is quite good even for 'less large' $a$: it has a maximum error of just $6\%$ near $a=3/2$. Of course for small $a$ the power series or Pade approximants are much better.

Answer 2

Using Mathematica:

p[x_, a_] = Map[a |-> PDF[NormalDistribution[a, 1], x], {-a, a}] // Mean // Simplify
h[a_] := NIntegrate[-p[x, a] Log[p[x, a]], {x, -Infinity, Infinity}]
v[x_, n_] := D[-p[x, a] Log[p[x, a]], {a, n}] /. a -> 0
c = Map[a |-> (x^a/a!) Integrate[v[x, a], {x, -Infinity, Infinity}], Range[0, 20, 2]] // Total
f[x_] = PadeApproximant[c, {x, 0, 4}]
Plot[{f[a], h[a]}, {a, -2, 2}]

I get some plot of the function between -2 and +2.

Here is the actual function: $$\frac{(1/2 (1 + \log(2 \pi)) + 1/2 x^2 (4 + 3 \log(2 \pi)) + 1/12 x^4 (22 + 7 \log(2 \pi)))}{(1 + 3 x^2 + (7 x^4)/6)}$$

Answer 3

Using your series for $h(a)$, let the constant apart.

What remains is an expansion starting with $a^2$; so the best Padé approximants will be $$h(a)=\frac{1+\log\left(2\pi\right)}{2}+\frac {a^2}2 \,\times\frac{1+\sum_{k=1}^n \alpha_k\,a^{2k}}{1+\sum_{k=1}^n \beta_k\,a^{2k}}$$

For example,

• using $n=2$, the difference with your series is $\color{red}{\large\frac{4933 }{25200}a^{12}}$

• using $n=3$, the difference with your series is $\color{red}{\large\frac{116680901 }{124311600}a^{16}}$

For $n=2$, the rational fraction is $$\frac {1+\frac{177 }{70}a^2+\frac{1}{210}a^4 } {1+\frac{106 }{35}a^2+\frac{83 }{70}a^4 }$$ and, for $n=3$, $$\frac {1+\frac{65027 }{9866}a^2+\frac{856039 }{103593}a^4-\frac{218397 }{690620}a^6 } {1+\frac{34980 }{4933}a^2+\frac{396266 }{34531}a^4+\frac{1799716 }{517965}a^6 }$$