# Computing Triple Integrals of Splines Fast

I'm trying to replicate a result from this paper, specifically in page 22, the function

$$I_1(s,s') = \max_{\phi \geq 2} \iiint\limits_{1/s' \leq t\leq u \leq v \leq 1/s} \omega \left( \frac{\phi - t- u - v}{u} \right) \frac{dt \, du \, dv}{tu^2v}$$

where $\omega(u)$ is the Buchstab function, defined by the following delay differential equation.

\begin{aligned} &\omega(u) = 1/u && 1 \leq u \leq 2 \\ &(u \omega(u))' = \omega(u-1) && u \geq 2 \end{aligned}

I already have the Buchstab function implemented as a piecewise set of Taylor polynomials $\omega_n(u)$, one for each interval $[n,n+1]$. Directly trying to compute the integration naively gets the computer stuck. Numerically integrating over the region in question seems to work:

NIntegrate[Boole[t <= u <= v] w[(p - t - u - v)/u]*1/(t*u^2*v),
{t, 1/sp, 1/s}, {u, 1/sp, 1/s}, {v, 1/sp, 1/s}]


but it's painfully slow. Any tips for getting this to run in reasonable time? Would numerically solving the differential equation and then trying to integrate that be better?

The definition of the Buchstab function I'm using follows from this paper (pages 137-138). We compute $\omega(u)$ by defining

$$\omega_j(u) = \omega(u), \; j \leq u \leq j+1$$

We have the following relation

$$u\omega_{j+1}(u) = \int_j^{u-1} \omega_j(t) dt + (j+1)\omega_j(j+1)$$

For each $\omega_j$, we expand it into a power series about $u=j+1$.

$$\omega_j(u) = \sum_{k=0}^\infty a_k(j) (u-(j+1))^k$$

We have $\omega_2(u)$ in exact form by the above relation, but the rest are hard to determine in closed form, as the repeat integrations get worse and worse.

$$\omega_2(u) = \frac{\log(u-1) + 1}{u}$$

I determined that the $a_k(j)$ coefficients are defined recursively,

\begin{align} a_0(2) &= \frac{1 + \log 2}{3} \\ a_k(2) &= (-1)^{k+1} \left( -\frac{1+\log 2}{3^{k+1}} + \frac{1}{3(2^k)} \sum_{m=0}^{k-1} \frac{1}{k-m} \left( \frac{2}{3}\right)^m \right) \\ a_0(j) &= a_0(j-1) + \frac{1}{j+1} \sum_{k=1}^\infty \frac{(-1)^k}{k+1} a_k(j-1) \\ a_k(j) &= \frac{a_{k-1}(j-1) - ka_{k-1}(j)}{k(j+1)} \end{align}

So I implement all this in Mathematica, and check it against the numerical solution to the DDE, and it all matches up.

wN[u_] := Evaluate[w[u] /.
NDSolve[{D[u*w[u], u] == w[u - 1], w[u /; u <= 2] == u^-1},
w, {u, 1, 30}, WorkingPrecision -> 50]];

w2Exact[u_] := (Log[u - 1] + 1)/u;

degreeTaylor = 10;

w2A[u_] := Evaluate[N[Normal[Series[w2Exact[u], {u,3,degreeTaylor}]],
MachinePrecision]];

w2Coeff = CoefficientList[w2A[u + 3], u];

a[k_, 2] := a[k, 2] = w2Coeff[[k + 1]];

a[0, j_] := a[0, j] = a[0, j - 1] + 1/(j + 1) Sum[(-1)^k/(k + 1)
*a[k,j - 1], {k, 1, degreeTaylor}];

a[k_, j_] := a[k, j] = 1/(j + 1) (a[k - 1, j - 1]/k - a[k - 1, j]);

w3A[u_] := Evaluate[Sum[a[k, 3]*(u - 4)^k, {k, 0, degreeTaylor}]];
w4A[u_] := Evaluate[Sum[a[k, 4]*(u - 5)^k, {k, 0, degreeTaylor}]];

wA[u_] := Piecewise[{
{u^-1, 1 <= u < 2},
{w2A[u], 2 <= u < 3},
{w3A[u], 3 <= u < 4},
{w4A[u], 4 <= u < 5},
{Exp[-EulerGamma], u >=  5}}];


Here wN is the numerical solution, and wA is my approximation. Plotting the two shows my approximation is good, and can be made as precise as needed.

LogPlot[Abs[wA[u] - wN[u]], {u, 1, 6}];


(I'm aware there is probably a better way to do it instead of having all these separate wnA functions, but trying something of the form wA[u_,n_] had issues I couldn't resolve.)

I add in as many $w_j$ terms as need be (for examples sake I'm only going up to $\omega_j(4)$ here), and increase the degree of the power series (here, only 10 terms) and then set it to its limiting constant $e^{-\gamma}$ after some point to obtain the accuracy needed.

It was done in this way so that the absolute error of my approximation to Buchstab could be computed in a rigorous manner, following the paper.

• For reference, could you include the definition you're using for the Buchstab function? And yes, NDSolve[] supports DDEs, so you might want to try that, too. Apr 4, 2017 at 8:26
• I'm having an awful lot of trouble formatting it, I've indented the code correctly and such, but I keep getting kicked back by the automated formatting checker, and it doesn't tell me where it's got a problem. Apr 4, 2017 at 9:37
• For now, can you try pasting it into Pastebin, and then somebody else can try formatting it into this post? Apr 4, 2017 at 9:39
• That's a good idea, pastebin.com/7S4WDc9r Apr 4, 2017 at 9:42
• P.S. Mathematics of Computation is freely accessible from AMS; I've replaced your JStor link because of this. Apr 4, 2017 at 9:51

The Buchstab function is the unique continuous function $w : \mathbb{R}_{\ge 1} \to \mathbb{R}_{>0}$ defined by the delay differential equation $u\,w(u) = 1 \; \, (1 \le u \le 2)$, $(u\,w(u))' = w(u - 1) \; \, (u > 2)$ (Panario, 1998). It approaches the asymptotic value $w(u) \to e^{-\gamma} \approx 0.561459$ as $u \to \infty$ (and in fact has nearly reached this value already by $u \approx 4$).

In The Mathematica GuideBook for Numerics | Michael Trott | Springer propose the following numerical solution, which I have taken care to copy in Mathematica:

BuchstabFunctionList =
FixedPointList[{(* integration interval start point *) #[[1]] + 1,
NDSolve[{D[z w[z], z] == #[[2]][[1, 1, 2]][z - 1],
(* solve differential equation *)
w[#[[1]] + 1] == #[[2]][[1, 1, 2]][#[[1]] + 1]}, w,
{z, #[[1]] + 1, #[[1]] + 2},
WorkingPrecision -> 40 , MaxSteps -> 10^4,
PrecisionGoal -> 20, AccuracyGoal -> 20]} &,
(* solution in interval [1, 2] *)
{1, {{w[z] -> (1/# &)} }},
(* solve until values at endpoints coincide *)
SameTest -> ((#1[[2, 1, 1, 2]][#1[[1]] + 1] ==
#2[[2, 1, 1, 2]][#2[[1]] + 1]) &)];

BuchstabFunction[x_] :=
BuchstabFunctionList[[IntegerPart[x]]][[2, 1, 1, 2]][x];

Plot[{BuchstabFunction[x], Exp[-EulerGamma]},
{x, 1, Length[BuchstabFunctionList] + 1},
AxesOrigin -> {1, 1},
AxesLabel -> {u, w[u]},
PlotRange -> {{1, 4}, {0.5, 1}},
PlotStyle -> {Red, {Blue, Dashed}}]


See if it can come in handy.

• Thankfully, solving DDEs is easier these days, compared to the time Trott wrote that code: NDSolve[{u ω'[u] + ω[u] == ω[u - 1], ω[u /; u <= 2] == 1/u}, ω, {u, 1, 30}]. Apr 4, 2017 at 8:58
• In this case, the DDE docs gives some nice examples. Apr 4, 2017 at 9:07
• The interval $[1,2]$ is implicitly taken care of by ω[u /; u <= 2] == 1/u and starting the integration at 1. Apr 4, 2017 at 9:19