I have had a similar question from my friend. He was surprised to find out $\displaystyle \lim_{x^{\circ} \rightarrow 0} \frac{\sin(x^{\circ})}{x^{\circ}} = \frac{\pi}{180}$.
The reason is that we are actually overloading the function name $\sin$. I think the confusion would be subsided if you were to look at $\sin_r(x)$ and $\sin_d(x)$ as two different functions where in the first function you input $x$ in radians and in the second function you input $x$ in degrees and these two are related by $\sin_r(x^r) = \sin_d(y^{\circ})$, where $y^{\circ} =\frac{180}{\pi}x^r$
We know that $\displaystyle \int_{0}^{\pi/2} \sin_r(x^r) dx^r = 1$.
Hence, $\displaystyle \int_{0}^{\pi/2} \sin_d(y^{\circ}) dx^r = 1$, since $\sin_r(x^r) = \sin_d(y^{\circ})$
We have $(\frac{180}{\pi}x^r) = y^{\circ}$. Hence $(\frac{180}{\pi}dx^r) = dy^{\circ}$.
As $x^r$ goes from $0$ to $\pi/2$, $y^{\circ}$ goes from $0$ to $90$.
Hence, we now have $\displaystyle \int_{0}^{90} \sin_d(y^{\circ}) \frac{\pi}{180} dy^{\circ} = 1$.
Hence, $\displaystyle \int_{0}^{90} \sin_d(y^{\circ}) dy^{\circ} = \frac{180}{\pi}$.
EDIT
The question is
"Why is $\int_{0}^{\pi/2} \sin(\theta) d\theta = 1$ when $\theta$ is in radians?"
What follows is my attempt for a purely geometric argument for this.
The geometric argument ultimately hinges (as one would expect) on
"The length of an arc of a circle subtending an angle $\theta$ at the center is $r \theta$ where $\theta$ is in radians"
This essentially comes from the way a radian is defined.
The main crux of the problem is the question
"Why is the average value of $\sin$ over quarter its period is $\frac{2}{\pi}$?"
(Note when we talk of average value of $\sin$ over quarter its period, there is no difference between radians or degrees).
It is easy to see that the average value of $\sin$ over quarter its period lies between $0$ and $1$.
The average value of of $\sin$ over quarter its period is nothing but the average value, with respect to $\theta$, of the height of the line segment as the line segment moves from the right end point of the circle towards the origin such that it is always perpendicular to the $X$ axis, with one end point of the segment on the $X$ axis and the other end point of the segment on the circumference of the circle.
Now here comes the claim. The average value, with respect to $\theta$, of the height of the line segment times the circumference of the quarter of the circle is the area of a unit square.
Let $h$ denote the average value, with respect to $\theta$, of the height of the line segment i.e. $$\displaystyle h = \frac{\int_{0}^{\pi/2} y d \theta}{\int_{0}^{\pi/2} d \theta}$$
The claim is
$$\displaystyle h \times \frac{\pi}{2} = 1$$
Proof:
First note that the area of the square with vertices at $(0,0),(0,1),(1,0),(1,1)$ is obviously $1$.
$\displaystyle \text{Area of square} = 1 = \int_{0}^{1} 1 \times dx = \int_{0}^{1} dx$.
For points on the circle, $\displaystyle x^2 + y^2 = 1 \Rightarrow xdx + ydy = 0 \Rightarrow dy = -\frac{x}{y} dx$.
Further, for the small triangle, we have
$$\displaystyle (dx)^2 + (dy)^2 = (1 \times d \theta)^2$$
$$\displaystyle (dx)^2 + \frac{x^2}{y^2}(dx)^2 = (d \theta)^2$$
$$\displaystyle \frac{x^2 + y^2}{y^2}(dx)^2 = (d \theta)^2$$
$$\displaystyle (dx)^2 = y^2 (d \theta)^2$$
Note that as $\theta$ increases from $0$ to $\frac{\pi}{2}$, $x$ decreases from $1$ to $0$ and hence
$$\displaystyle dx = -y d \theta$$
Hence, $$\displaystyle \text{Area of square} = 1 = \int_{0}^{1} 1 \times dx = \int_{0}^{1} dx = \int_{\pi/2}^0 (-y) d \theta = \int_0^{\pi/2} y d \theta$$
Hence, $\displaystyle 1 = \int_0^{\pi/2} y d \theta = \frac{\int_0^{\pi/2} y d \theta}{\int_0^{\pi/2} d \theta} \times \int_0^{\pi/2} d \theta$
Note that $\displaystyle \int_0^{\pi/2} d \theta = \frac{\pi}{2}$ which is the circumference of the quadrant and $\displaystyle h = \frac{\int_0^{\pi/2} y d \theta}{\int_0^{\pi/2} d \theta}$ is the average value, with respect to $\theta$, of the height of the line segment.
Hence, we get $\displaystyle h \times \frac{\pi}{2} = 1$ i.e. average value, with respect to $\theta$, of the height of the line segment if $\frac{2}{\pi}$.
Hence, we get that $$\displaystyle \int_{0}^{\pi/2} \sin(\theta) d \theta = 1$$ when $\theta$ is in radians
I tried 2 ways to find a closed form, though unsuccessful up to this point.
1st trial. Let $I$ denote the integral, and write
$$ I = 8 \sum_{n=0}^{\infty} \frac{(-1)^{n}}{2n+1} \int_{0}^{\infty} \frac{e^{-x}(1 - e^{-(2n+1)x})}{x (1 + e^{-x})^{2}} \, dx. $$
In order to evaluate the integral inside the summation, we introduce new functions $I(s)$ and $J_n(s)$ by
$$ I(s) = 8 \sum_{n=0}^{\infty} \frac{(-1)^{n}}{2n+1} \int_{0}^{\infty} \frac{x^{s-1} e^{-x}(1 - e^{-(2n+1)x})}{(1 + e^{-x})^{2}} \, dx =: 8 \sum_{n=0}^{\infty} \frac{(-1)^{n}}{2n+1} J_n(s) $$
so that $I = I(0)$. Then it is easy to calculate that for $\Re(s) > 1$, $J_n(s)$ is written as
$$ J_n(s) = \Gamma(s) \left( \eta(s-1) + \sum_{k=2n+1}^{\infty} \frac{(-1)^{k-1}}{k^{s-1}} - (2n+1) \sum_{k=2n+1}^{\infty} \frac{(-1)^{k-1}}{k^{s}} \right), $$
where $\eta$ is the Dirichlet eta function. Plugging this back to $I(s)$ and manipulating a little bit, we obtain
$$ I(s) = 8\Gamma(s) \left( \frac{\pi}{4} \eta(s-1) - 4^{-s}\left( \zeta(s, \tfrac{1}{4}) - \zeta(s, \tfrac{1}{2}) \right) + \sum_{n=0}^{\infty} \frac{(-1)^{n}}{2n+1} \sum_{k=2n+1}^{\infty} \frac{(-1)^{k-1}}{k^{s-1}} \right). $$
This is valid for $\Re(s) > 1$. But if we can somehow manage to find an analytic continuation of the last summation part, then we may find the value of $I = I(0)$.
2nd trial. I began with the following representation
\begin{align*}
I
&= -2 \int_{0}^{\infty} \frac{1-e^{-t}}{1+e^{-t}} \left( \frac{1}{\cosh t} - \frac{2}{t} ( \arctan(1) - \arctan (e^{-t})) \right) \, \frac{dt}{t} \\
&= \sum_{n=0}^{\infty} \frac{(-1)^{n}}{(2n+1)(2n+2)} \int_{-\infty}^{\infty} \frac{\tanh^{2(n+1)} x}{x^{2}} \, dx.
\end{align*}
With some residue calculation, we can find that
\begin{align*}
\int_{-\infty}^{\infty} \frac{\tanh^{2n} x}{x^{2}} \, dx
&= \frac{2}{i\pi} \, \underset{z=0}{\mathrm{Res}} \left[ \psi_{1}\left(\tfrac{1}{2} + \tfrac{1}{i\pi} z\right) \coth^{2n} z \right] \\
&= 2^{2n+3} \sum_{m=1}^{n} (-1)^{m-1}m (1-2^{-2m-1}) A_{n-m}^{(2n)} \, \frac{\zeta(2m+1)}{\pi^{2m}},
\end{align*}
where $A_m^{(n)}$ is defined by the following combinatoric sum
$$ A_m^{(n)} = \sum_{\substack{ k_1 + \cdots + k_n = m \\ k_1, \cdots, k_n \geq 0 }} \frac{B_{2k_1} \cdots B_{2k_n}}{(2k_1)! \cdots (2k_n)!} = 2^{-2m} [z^{2m}](z \coth z)^{n} \in \Bbb{Q}, $$
where $B_k$ are Bernoulli numbers. Still the final output is egregiously complicated, so I stopped here.
3rd trial. The following yet another representation may be helpful, I guess.
$$ I = \int_{0}^{1/2} \frac{1 - \cot(\pi u/2)}{2} \left\{ \psi_1\left(\tfrac{1+u}{2}\right) - \psi_1\left(\tfrac{1-u}{2}\right) \right\} \, du. $$
Best Answer
Answering a comment by the OP :
You can improve the upper bound for example by finding the tangent of $f(x)=-x^2$ with slope $-1$. Then $f'(x_0) = -2x_0 =-1$ , so it is the tangent at $x_0 = \frac{1}{2}$ , its function is $y + \frac{1}{4} = -(x-\frac{1}{2}) $ which gives you $y = -x + \frac{1}{4}$ which means that $\int_0^1 e^{-x^2} < \int_0^1 e^{-x+\frac{1}{4}}= 0.811659... $;
You can also use a taylor approximation to find better estimates indeed the function $e^{-x^2}$ is bounded below by a taylor approximation with an even number of terms and above by a taylor approximation with an odd number of terms
$$ 1-x^2+\frac{x^4}{2!}-...+(-1)^{2k+1}\frac{x^{2(2k+1)}}{(2k+1)!}< e^{-x^2} < 1-x^2+\frac{x^4}{2!}-...+ (-1)^{2k} \frac{x^{2(2k)}}{(2k)!} $$