# Integrating exponential of second degree polynomials

I'm looking to compute the value of the following integral, for small values of $$|a|$$. $$u_n(a,b)=\frac{1}{2}\int_{-1}^1 x^ne^{ax^2+bx}\mathrm{d}x$$

In this equation, $$a,b \in \mathbb{R}$$ and $$n \in \mathbb{N}$$.

It's easy to show that:

• $$u_0(0, 0) = 1$$
• $$\frac{\partial}{\partial b} u_n = u_{n+1}$$
• $$\frac{\partial^2}{\partial b^2} u_n = \frac{\partial}{\partial a} u_n = u_{n+2}$$

We can find analytic solutions for this equation using Mathematica, for example here is $$u_0$$: $$u_0(a,b) = e^a \frac{e^b F\left(\frac{b+2a}{2\sqrt a}\right)-e^{-b} F\left(\frac{b-2a}{2\sqrt a}\right)}{2\sqrt a}$$ where $$F(x)$$ is the Dawson function. We can also show that: $$\lim_{a\to0} u_0 = \frac{\sinh(b)}{b}$$ $$\lim_{b\to0} u_0 = e^a\frac{F(\sqrt a)}{\sqrt a}$$

The analytic solution works great for larger $$|a|$$ as there are a few libraries that can provide the value of the Dawson function. However, for small values of $$|a|$$, the analytic solution has a catastrophic $$\frac{0}{0}$$ limit. This problem only gets worse as we go to $$u_1$$, $$u_2$$, etc.

Is there a way to reformulate this problem in a way where it can be computed for small values of $$|a|$$? For example, perhaps the value of $$u_n(a, b)$$ be the limit of some iteration?

The best I could do for now is to compute a Taylor series around $$a=0$$: $$u_n(a,b) = \sum_k \frac{a^k}{k!} v_{2k+n}(b)$$ where $$v_0 = \frac{\sinh(b)}{b}$$ and we use the recurrence relation: $$b v_{n+1} = f(b) - (n+1) v_n$$ Here $$f(b)$$ is $$\cosh(b)$$ if $$n$$ is even and it is $$\sinh(b)$$ if $$n$$ is odd. However, to compute for $$a = 0.5$$, we need quite a few terms and it's slow to compute.

Edit:

Here is a small induction proof for the recurrence relation. We simply derive it by $$b$$ and note that $$\frac{\partial}{\partial b}v_n = v_{n+1}$$ as $$v_n = \lim_{a\to0}u_n$$ and $$\frac{\partial}{\partial b}u_n = u_{n+1}$$. For $$n$$ even, we have: $$\frac{\partial}{\partial b} bv_{n+1} = \frac{\partial}{\partial b} \left( \cosh(b) - (n+1)v_n\right)$$ $$bv_{n+2} + v_{n+1} = \sinh(b) - (n+1)v_{n+1}$$ $$bv_{n+2} = \sinh(b) - (n+2)v_{n+1}$$

We have a similar proof for the case $$n$$ is odd (simply exchange $$\cosh$$ and $$\sinh$$).

• So, for a=0.5, why not use that analytic formula with the Dawson function? Nothing catastrophic should be happening there at a=0.5. Sep 28, 2021 at 14:52
• The problem arises for larger values of $n$. For example, we can show that: $$u_1 = \frac{b}{2a}\left(e^a \frac{\sinh(b)}{b} - u_0\right)$$ In this case, we have one more order of the problematic $\frac{0}{0}$ that appears as $a \to 0$.
– PC1
Sep 28, 2021 at 15:19

By inspecting more carefully nicoguaro's solution, I realize that in general: $$v_n(b) = (-1)^{n+1} \frac{\Gamma(n+1,b) - \Gamma(n+1,-b)}{2b^{n+1}}$$

It's straightforward to confirm that $$v_0 = \frac{\sinh(b)}b$$ and that $$\frac{\partial}{\partial b} v_n = v_{n+1}$$ for all $$n \in \mathbb{N}$$, so this is the same function!

So we can use the Taylor expansion (to any arbitrary order in $$a$$) with the recurrence relation as we can directly compute any starting $$v_n$$ term and then propagate to nearby values of $$v_{n+2}$$, $$v_{n+4}$$, etc. using the recurrence relation provided in the original question.

• @nicoguaro FYI I think that this solution can work, thank you for your input!
– PC1
Sep 29, 2021 at 16:13

I think that a numerical integration might work in your case, I don't see why it should not.

Besides numerically integrating your function, you could try an asymptotic expansion for small $$a$$. Since you are using Mathematica, I tried to compute the integral using Wolfram Language with the function AsymptoticIntegrate, specifically I used

AsymptoticIntegrate[x^n * Exp[a*x^2 + b*x], {x, -1, 1}, a->0]


This gave me the following result that seems to be fine for the values I tested.

$$u_n(a, b) \sim \frac{1}{b^3}\left(-b^2\right)^{-n} \left[b^2 \left((-1)^n (-b)^n-b^n\right) \Gamma(1+n)+a \left((-1)^n (-b)^n-b^n\right) \Gamma(3+n)+b^{2+n} \Gamma(1+n,-b)+(-1)^{1+n} (-b)^{2+n} \Gamma(1+n,b)+a b^n \Gamma(3+n,-b)+(-1)^{1+n} a (-b)^n \Gamma(3+n,b)\right]$$

• Thank you for the suggestion. I am not too familiar with the incomplete gamma function, are you aware of any software (for example in C or C++) that can efficiently compute it for large values of $n$? Typically, I will need to compute it for $20 < n < 50$.
– PC1
Sep 28, 2021 at 21:05
• @PC1, according to GSL they can handle it up to 171. Sep 29, 2021 at 10:52