How can you calculate $\int_{0}^{1} \exp(\sin(x))dx$ exactly?

I'm taking a computational statistics class and I have to calculate $\int_{0}^{1} \exp(\sin(x))dx$ numerically (Monte Carlo integration) and compare the numerical result to the exact result. The problem is, how can one canculate the exact value of the integral? The antiderivative doesn't seem to exist in terms of elementary functions. Mathematica gives me this:

Not only that, Mathematica also cannot calculate the definite integral:

And since Mathematica cannot do it, I don't think other softwares can, though I'm not sure. Of course, Mathematica can do the numerical integration

But that's not what I need! Thanks in advance for any help.

To integrate

$$\int_0^1 e^{\sin(x)}\ dx$$

take first $x = \sin^{-1}(u)$, so $dx = \frac{1}{\sqrt{1 - u^2}}\ du$, then

$$\int_{0}^{1} e^{\sin(x)}\ dx = \int_{0}^{\sin(1)} \frac{e^u}{\sqrt{1 - u^2}}\ du$$

and now note that, by using the general binomial theorem and the identity $(m)_k = (-1)^k (-m)^{(k)}$ converting the falling and rising factorial power,

$$\begin{align}\frac{1}{\sqrt{1 - u^2}} &= (1 - u^2)^{-1/2} \\ &= \sum_{k=0}^{\infty} \binom{-1/2}{k} (-1)^k u^{2k} \\ &= \sum_{k=0}^{\infty} \frac{(-1/2)_k}{k!} (-1)^k u^{2k}\\ &= \sum_{k=0}^{\infty} (1/2)^{(k)} \frac{(u^2)^k}{k!}\end{align}$$

so that, by simply multiplying the series, we have

$$\frac{e^u}{\sqrt{1 - u^2}} = \sum_{n=0}^{\infty} \sum_{k=0}^{\infty} (1/2)^{(k)} \frac{u^n}{n!} \frac{(u^2)^k}{k!}$$

And integrating, we get, noting that we should combine $u^n$ and $(u^2)^k = u^{2k}$ to $u^{(2k + n)}$,

$$\int_{0}^{\sin(1)} \frac{e^u}{\sqrt{1 - u^2}}\ du = \int_{0}^{1} e^{\sin(x)} dx = \sum_{n=0}^{\infty} \sum_{k=0}^{\infty} \frac{(1/2)^{(k)}}{(2k + n + 1)} \frac{[\sin(1)]^{2k + n + 1}}{n! k!}$$

More generally,

$$\int e^{\sin(x)}\ dx = C + \sum_{n=0}^{\infty} \sum_{k=0}^{\infty} \frac{(1/2)^{(k)}}{(2k + n + 1)} \frac{[\sin(x)]^{2k + n + 1}}{n! k!}$$

. Note that this is a hypergeometric series in two variables. To see that, you should first note the identity $\frac{2^{(n)}}{1^{(n)}} = n + 1$, thus

$$\sum_{n=0}^{\infty} \sum_{k=0}^{\infty} \frac{(1/2)^{(k)}}{(2k + n + 1)} \frac{u^{2k + n + 1}}{n! k!} = u \sum_{n=0}^{\infty} \sum_{k=0}^{\infty} \frac{(1/2)^{(k)} 1^{(2k+n)}}{2^{(2k+n)}} \frac{(u^2)^k u^n}{k!n!}$$

so it looks like some sort of

$$F_\mathrm{hyp}(a; b; c; x, y) := \sum_{n=0}^{\infty} \sum_{k=0}^{\infty} \frac{a^{(n)} b^{(2n+k)}}{c^{(2n+k)}} \frac{x^n y^k}{n! k!}$$

double hypergeometric function in $u = \sin(x)$ and $[\sin(x)]^2$. But I don't think this particular form has any official name, though with some work maybe it can be reduced to known ones. The Horn functions 1 come closest, with the doubled index; but none have a $2n+k$ index in both numerator and denominator of the coefficient part. Two-variable hypergeo are obscure, though, at least from Wolfram's point of view, because it does not implement them.

Welcome to this new and exotic world.

Since $1$ is not a natural bound of integration for this problem, it is pretty clear that you would have to have an indefinite integral for this to have any chance of success. By contrast, if the upper limit were $\pi,$ you get the nice answer:

$$\int_{0}^{\pi} \exp(\sin(x))dx = \pi (\pmb{L}_0(1)+I_0(1)),$$ while for $2\pi$ you get an even nicer answer:

$$\int_{0}^{2\pi} \exp(\sin(x))dx = 2 \pi I_0(1).$$

(with $I$ being the Bessel function and $L$ the Struve $L$ function.

As @Bumblebee commented, for the definite integral, compose Taylor series for the integrand $$e^{\sin(x)}=1+x+\frac{x^2}{2}-\frac{x^4}{8}-\frac{x^5}{15}-\frac{x^6}{240}+\frac{x^7}{90}+\frac{31 x^8}{5760}+\frac{x^9}{5670}-\frac{2951 x^{10}}{3628800}-\frac{x^{11}}{3150}+\frac{181 x^{12}}{14515200}+\frac{2417 x^{13}}{48648600}+\frac{58913 x^{14}}{4151347200}-\frac{5699 x^{15}}{2554051500}+O\left(x^{16}\right)$$

Integrate termwise and use the bounds to get $$\int_0^1 e^{\sin(x)}\,dx\sim \frac{60970121917}{37362124800}=1.63186977\cdots$$ while numerical integration gives $1.63186961$

Edit

The next term in the expansion is $-\frac{52635599 }{20922789888000}x^{16}$ that is to say that its contribution to the integral is $-\frac{52635599}{355687428096000} \sim -1.48\times 10^{-7}$ which is not much.

Update (almost eight years later)

Using $$e^{\sin(x)}=\sum_{n=0}^\infty \frac{\sin^n(x)}{n!}$$ $$I_n=\int \sin^n(x)\,dx=-\cos (x) \, _2F_1\left(\frac{1}{2},\frac{1-n}{2};\frac{3}{2};\cos ^2(x)\right)$$ where appears the gaussian hypergeomtric function. $$J_n=\int_0^t \sin^n(x)\,dx=\frac{\sqrt{\pi }\,\, \Gamma \left(\frac{n+1}{2}\right)}{2 \Gamma \left(\frac{n+2}{2}\right)}-\cos (t) \, _2F_1\left(\frac{1}{2},\frac{1-n}{2};\frac{3}{2};\cos ^2(t)\right)$$

For $t=1$, the partial sums converge very fast $$\left( \begin{array}{cc} p & \sum_{n=0}^p \frac {J_n}{n!} \\ 0 & 1.0000000 \\ 1 & 1.4596977 \\ 2 & 1.5960355 \\ 3 & 1.6258589 \\ 4 & 1.6310267 \\ 5 & 1.6317681 \\ 6 & 1.6318589 \\ 7 & 1.6318686 \\ 8 & 1.6318695 \\ 9 & 1.6318696 \end{array} \right)$$