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.


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$).

  • $\begingroup$ 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. $\endgroup$ Sep 28, 2021 at 14:52
  • 1
    $\begingroup$ 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$. $\endgroup$
    – PC1
    Sep 28, 2021 at 15:19

2 Answers 2


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.

  • $\begingroup$ @nicoguaro FYI I think that this solution can work, thank you for your input! $\endgroup$
    – 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]$$

  • $\begingroup$ 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$. $\endgroup$
    – PC1
    Sep 28, 2021 at 21:05
  • 1
    $\begingroup$ @PC1, according to GSL they can handle it up to 171. $\endgroup$
    – nicoguaro
    Sep 29, 2021 at 10:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.