# 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. – J. M. will be back soon Apr 4 '17 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. – J. Ashford Apr 4 '17 at 9:37
• For now, can you try pasting it into Pastebin, and then somebody else can try formatting it into this post? – J. M. will be back soon Apr 4 '17 at 9:39
• That's a good idea, pastebin.com/7S4WDc9r – J. Ashford Apr 4 '17 at 9:42
• P.S. Mathematics of Computation is freely accessible from AMS; I've replaced your JStor link because of this. – J. M. will be back soon Apr 4 '17 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}]. – J. M. will be back soon Apr 4 '17 at 8:58
• @ J.M.: Before I die I have to find out how do you always find the right code. I have tried repeatedly but unlike your writing 1 <= u <= 2 and so I did not get anything. : '( This is why I decided to copy that code because I do not think he did it anyone before. – TeM Apr 4 '17 at 9:02
• In this case, the DDE docs gives some nice examples. – J. M. will be back soon Apr 4 '17 at 9:07
• I have referred precisely to what is written there in my attempts. But I can not understand why you can not see the interval closed and limited [1, 2]. – TeM Apr 4 '17 at 9:10
• The interval $[1,2]$ is implicitly taken care of by ω[u /; u <= 2] == 1/u and starting the integration at 1. – J. M. will be back soon Apr 4 '17 at 9:19