It seems, the main difficulty in defining integration on surreals stems from the set of axioms the authors decided to choose for their definition of integration to satisfy.
Particularly, the most problematic is the property of linearity against an infinitely-large factor. For instance, it is the property $(b)$ in Preposition 14 in this paper on surreal integration:
$$\int_a^b (\alpha f+\beta g)=\alpha \int_a^b f + \beta \int_a^b g,$$
where $\alpha, \beta \in \text{No}$.
In my view, this property is unnatural. For instance, assuming total order, you definitely want a surreal number, such big that it to correspond to the derivative of a function that makes a jump at a point, like $\operatorname{sgn}(x)$ makes at zero. You definitely want a surreal number such big that integration of a function over a range where the function takes this number as a value only at one point is non-zero.
Otherwise, you would be able to extend your order of surreal numbers by new numbers who would also be totally ordered. But surreal numbers are a maximal ordered field.
Once this axiom is dropped and replaced by the linearity against finite (real) values, the integration can be easily defined based on this natural definition:
$$\int_S u dt=\pi N(S) \partial(u),$$
where $u$ is a surreal constant, $N(S)$ is numerosity and $\partial(u)$ is derivation. The numerosity plays here the same role as measure in Lebesgue integral. Using this base formula, one can derive a formula for integrating arbitrary surreal-valued functions:
$$\int_a^b F(t) dt=\pi \omega_1\tilde{\int}_a^b \partial{(F(t))}dt+\int_a^b \operatorname{fin}F(t)dt,$$
where $b>a; a,b\in\text{No}$, $F$ is a surreal-valued function of surreal argument, the symbol $\tilde{\int}$ should mean that the integral is taken formally, without awareness of the surreal character of the expression under the integral, and $\operatorname{fin}F(t)$ means taking finite part of $F(t)$. The second integral is evaluated using Newto-Leibnitz formula after finding the indefinite integral symbolically.
Below is the code for Mathematica that realizes this formula. For simplicity, in input w
corresponds to $\omega$ and W
to $\omega_1$. The input accepts any surreal numbers composed of these values as integration limits and any surreal-valued function of integration variable x
as the expression under integral. In output, derivation is denoted as ꝺ
.
f[x_] = w^W;
a = 0;
b = w;
f[x_] = f[x] /. W -> W[w];
Unprotect[Power]; Power[0, 0] = 1; Protect[Power];
Fin[t_] := (PowerExpand[f[t]] /. Log[w] -> 0 /. w -> 0) +
Limit[Evaluate[LaplaceTransform[ D[f[t], w], t, x]] +
Evaluate[LaplaceTransform[ D[f[t], w], t, -x]], x -> 0]/
2 /. \[Infinity] -> 0 /. -\[Infinity] -> 0 //
FullSimplify (*Finding the finite part so to integrate separately*)
int = \[Pi] W Integrate[D[f[t], w], {t, a, b}] +
Integrate[Fin[t], {t, a, b}] // FullSimplify;
int = If[a == b, \[Pi] D[f[t], w], int, int] /.
Derivative[1][W][w] -> ꝺ[W] /. W[w_] :> W // Normal;
Print[Inactivate[
Integrate[
f[x] /. W[w] -> W /. w -> \[Omega] /.
W -> Subscript[\[Omega], 1], {x,
a /. w -> \[Omega] /. W -> Subscript[\[Omega], 1],
b /. w -> \[Omega] /. W -> Subscript[\[Omega], 1]}],
Integrate], "=",
int /. w -> \[Omega] /. W -> Subscript[\[Omega], 1]]
The code correctly gives all the results from the other my answer, but also can find more complicated results, such as:
$\int _0^{\omega }e^{c x \omega }dx=\frac{c \omega ^3+\omega +\pi \left(e^{c \omega ^2} \left(c \omega ^2-1\right)+1\right) \omega _1}{c \omega ^2}$
$\int _0^{\omega }e^{c x^2 \omega }dx=\omega +\frac{1}{2} e^{c \omega ^3} \pi \omega _1-\frac{\pi ^{3/2} \text{erfi}\left(\sqrt{c} \omega ^{3/2}\right) \omega _1}{4 \sqrt{c} \omega ^{3/2}}$
$\int _0^{\omega }\log (c x \omega )dx=\omega (\log (c)+\log (\omega )-1)+\pi \omega _1$
$\int _0^{\omega }\omega ^xdx=\frac{\frac{\pi \left((\omega \log (\omega )-1) \omega ^{\omega }+1\right) \omega _1}{\omega }+1}{\log ^2(\omega )}$
$\int _0^{\omega }\omega ^{\omega }dx=\pi (\log (\omega )+1) \omega _1 \omega ^{\omega +1}+\omega$
$\int _0^{\omega }p \omega ^{q \omega }dx=p \omega \left(\pi q (\log (\omega )+1) \omega _1 \omega ^{q \omega }+1\right)$
$\int _0^{\omega }\omega ^{\omega _1}dx=\pi \omega _1 \left(\omega \partial \left(\omega _1\right) \log (\omega )+\omega _1\right) \omega ^{\omega _1}$
$\int _0^{\omega }\omega _1^{\omega }dx=\pi \omega \left(\frac{\omega \partial\left(\omega _1\right)}{\omega _1}+\log \left(\omega _1\right)\right) \omega _1^{\omega +1}+\omega$
$\int _0^{\omega }\omega _1^xdx=\frac{\omega _1^{\omega }-1}{\log \left(\omega _1\right)}+\frac{\partial\left(\omega _1\right) \left(\pi \left(\omega \log \left(\omega _1\right)-1\right) \omega _1^{\omega +1}+\pi \omega _1+\omega \right)}{\log ^2\left(\omega _1\right) \omega _1}$
$\int _0^{\omega _1}\omega ^xdx=\frac{\omega _1 \left(\pi \left(\log (\omega ) \omega _1-1\right) \omega ^{\omega _1}+\pi +1\right)}{\omega \log ^2(\omega )}$
$\int _0^{\omega }\log \left(\omega _1\right)dx=\omega \left(\log \left(\omega _1\right)+\pi \partial\left(\omega _1\right) \omega _1\right)$
$\int _0^{\omega }\log \left(\omega _1\right)dx=\omega \left(\log \left(\omega _1\right)+\frac{\pi \partial\left(\omega _1\right) \omega _1}{\omega _1(\omega )}\right)$
etc.