# How to speed up the integral in NDSolve?

On the MMA.SE, there have been several question on how to solve equations including integral using NDSolve, see example 1, example 2, and example 3. Many people, like me, have a hard time with this kind of problem. I had raised the last one 5 months ago and obtained a promising answer, since then I have been tried to solve the following equation for a spatially periodic function $$u(x,t)$$ on $$[-L,L]$$ with $$2L$$ periodicity:

$$\partial_t u + u\partial_x u + \partial_x^2 u +\partial_x^4 u + a{\rm{int}}[\partial_x^3u] + bu^3 {\rm{int}} [(\partial_x^2u) {\rm{int}}[\partial_x u] ] = 0, \tag{1}$$

where $$a$$ and $$b$$ are constants, and $$\rm{int}[f]$$ is a spatial integral for a periodic function $$f(x,t)$$

$${\rm{int}}[f](x,t)=\frac{1}{2L}\mathrm{PV}\int_{-L}^L f(x^\prime,t)\cot\left[\frac{\pi(x-x^\prime)}{2L}\right]\:\mathrm{d}x^\prime, \tag{2}$$

which should be understood in the sense of principal value (PV) since there is a singularity at $$x=x^\prime$$.

This equation is subjected to periodic boundary conditions and an initial condition. In my real problem, it has several $$\rm{int}$$ terms and its nest, like $${\rm{int}}[(\partial_x^2u) \rm{int}[\cdots] ]$$, which will be solved over a large domain and a long time, say, $$L=30$$ and $$t_\text{max}=200$$ (much larger than those in my original post). I tried Michael's answer but found that even for that simplified version of Eq.(1) (without the last nested term) the code is extremely slow, though it works. Actually, I also posted my answer based on the finite difference method, but the accuracy of my code is not as good as Michael's. (That is why I didn't accept any answer.) Also, the problem has stiffness due to the high-order derivatives and nonlinear terms. So I have added some Method options to NDSolve ProcessEquations.

The Mathematica code

L = 30; tmax = 30; a = 1; b = 1/100; c = 1/(2 L); e = 1/10; nGrid = 91;
ic[x_] = e*Cos[\[Pi]*x/L];

sys = {D[u[x, t], t] + u[x, t]*D[u[x, t], x] + D[u[x, t], {x, 2}] +
D[u[x, t], {x, 4}] + a*int[D[u[x, t], {x, 3}], x, t] +
b*u[x, t]^3*intnest[D[u[x, t], {x, 2}]*int[D[u[x, t], x], x, t], x, t] == 0, u[-L, t] == u[L, t], u[x, 0] == ic[x]};

periodize[data_] := Append[data, {N@L, data[[1, 2]]}];(*for periodic interpolation*)

Block[{int, intnest},
(* IC fools ProcessEquations to consider int[] as a good num.fn.*)
int[uppp_, x_?NumericQ, t_ /; t == 0] := (cnt++;
c*NIntegrate[D[ic[xp], {xp, 3}]*Cot[\[Pi] (x - xp)/(2*L)], {xp, x - L, x, x + L},
Method -> {"InterpolationPointsSubdivision", Method -> {"PrincipalValue", "SymbolicProcessing" -> 0}},
PrecisionGoal -> 8, AccuracyGoal -> 8, MaxRecursion -> 10]);
int[uppp_?VectorQ, xv_?VectorQ, t_?NumericQ] := Function[x, cnt++;
c*NIntegrate[Interpolation[periodize@Transpose@{xv, uppp}, xp,
PeriodicInterpolation -> True]*Cot[\[Pi] (x - xp)/(2*L)], {xp, x - L, x, x + L},
Method -> {"InterpolationPointsSubdivision", Method -> {"PrincipalValue", "SymbolicProcessing" -> 0}},
PrecisionGoal -> 8, AccuracyGoal -> 8, MaxRecursion -> 10]] /@xv;
intnest[upp_, x_?NumericQ, t_ /; t == 0] := (cnt2++;
c*NIntegrate[D[ic[xp], {xp, 2}]*int[D[ic[xp], xp], x, t]*Cot[\[Pi] (x - xp)/(2*L)], {xp, x - L, x, x + L},
Method -> {"InterpolationPointsSubdivision", Method -> {"PrincipalValue", "SymbolicProcessing" -> 0}},
PrecisionGoal -> 8, AccuracyGoal -> 8, MaxRecursion -> 10]);
intnest[upp_?VectorQ, xv_?VectorQ, t_?NumericQ] := Function[x, cnt2++;
c*NIntegrate[Interpolation[periodize@Transpose@{xv, upp}, xp,
PeriodicInterpolation -> True]*Cot[\[Pi] (x - xp)/(2*L)], {xp, x - L, x, x + L},
Method -> {"InterpolationPointsSubdivision", Method -> {"PrincipalValue", "SymbolicProcessing" -> 0}},
PrecisionGoal -> 8, AccuracyGoal -> 8, MaxRecursion -> 10]] /@xv;
(*monitor while integrating pde*)
Clear[foo];
cnt = 0; cnt2 = 0;
PrintTemporary@Dynamic@{foo, cnt, cnt2, Clock[Infinity]};
(*broken down NDSolve call*)
InternalInheritedBlock[{MapThread},
{state} = NDSolveProcessEquations[sys, u, {x, -L, L}, {t, 0, tmax},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> nGrid, "MaxPoints" -> nGrid, "DifferenceOrder" -> "Pseudospectral"},
Method -> {"StiffnessSwitching", "NonstiffTest" -> Automatic}},
AccuracyGoal -> Infinity, WorkingPrecision -> 20,
MaxSteps -> \[Infinity], StepMonitor :> (foo = t)];
Unprotect[MapThread];
MapThread[f_, data_, 1] /; ! FreeQ[f, int] := f @@ data;
Protect[MapThread];
NDSolveIterate[state, {0, tmax}];
sol = NDSolveProcessSolutions[state]]] // AbsoluteTiming


My problem

As mentioned above, the code is very slow. An estimation: $$>2$$ hrs may be required to obtain convergence with tmax = 1. Btw, the "slwcon" and "ncvb" warning could be ignored (see Michael's comments following his answer there). Is there any approach that would help speed up the code? Thank you very much.

Some ideas

As suggested by Henrik Schumacher, the combination of NIntegrate and Interpolation limits the speed of this code. Maybe it could be better to use a fixed quadrature rule and implement the integration with ListConvolve. But I need help with implementing this idea in my problem, so I bring this problem here in the hope that someone could help.

I am thinking that can we divide the interval into a uniform grid with $$2M$$($$=\rm{nGrid}-1$$) mesh points defined by $$x_m=(m-M)h$$, where $$h=L/M$$. Please see also my answer to a similar problem. Then the integral term (2) could be evaluate at the midpoints $$x_{i+1/2}=(x_i+x_{i+1})/2$$, for $$i=0,1,\ldots,2M-1$$, (note the periodicity demands $$u_0=u_{2M}$$) using a certain integration rule, e.g., trapezoidal rule, with $$x_i$$ as integration nodes. In this way, the principal value integral could be efficiently computed, as if it were simply an ordinary integral.

• This PDE is parabolic of order 5. Those are hard (and slow) to integrate, even with implicit methods. With the integral term: Mathematica does not know how to handle it implicitly. Moreover, the combination of NIntegrate and Interpolation is bound to be very slow. Maybe it s better to find a suitable fixed quadrature rule and implement it by hand and compute the integral by ListConvolve. Jul 17 '19 at 12:20
• @HenrikSchumacher could you show an example by computing the integral by ListConvolve using, say, trapezoidal integration or provide a relevant link? Thank you very much! Jul 18 '19 at 4:40
• @HenrikSchumacher Is it possible to compile the code (partially) in order to speed up the computation? Jul 18 '19 at 11:19
• No, integration of Integration does not help. Jul 18 '19 at 12:04
• Cimpilation of Integrate is useluess in_any_ case unless the integral can be computed analytically. Jul 19 '19 at 2:39

## 1 Answer

It is possible to solve this problem using the decomposition of the solution in a Fourier series. Then it is possible to replace the integrals with the coefficients of the Fourier series using the following obvious property
$$\frac {1}{2\pi}\int _{-\pi}^{\pi}e^{-inx'}\cot((x-x')/2)dx'=ie^{-inx}(1+ \mbox{sign}(n))$$ It is exactly equal to FourierCoefficient[Cot[(x - xp)/2], xp, n] The result is a fairly simple code. But the calculation does not coincide with that obtained using the author's code. It means that we have to write another third code to test these two. For ease of use of Fourier series, we convert the coordinate and time according to x->k0 x, t->k0 t, k0=Pi/L.

L = 30; tmax = 30; a = 1; b = 1/100; c = 1/(2 L); e = 1/10; nn = 10; k0 = Pi/L; tm = tmax*k0;
a1 = a/k0;
uf[x_, t_] := Sum[f[k][t] Exp[I k x], {k, -nn, nn}]
eq = Table[
f[m]'[t] +
I Sum[ If[Abs[m - k] <= nn, f[m - k][t], 0] k f[k][t], {k, -nn,
nn, 1}] - k0*m^2 f[m][t] + k0^3 m^4 f[m][t] -
a k0^2 (I m)^3 I (1 - Sign[m]) f[m][t] +
b k0^2 Sum[
If[Abs[m - k - s1 - s2 - s3] <= nn,
f[s1][t] f[s2][t] f[s3][t] f[k][t] f[m - k - s1 - s2 - s3][
t] I (1 - Sign[k]) I (1 - Sign[m - k - s1 - s2 - s3]),
0], {s1, -nn, nn}, {s2, -nn, nn}, {s3, -nn, nn}, {k, -nn,
nn}] == 0, {m, -nn, nn, 1}];

ic = Table[
f[m][0] ==
e (KroneckerDelta[m, 1] + KroneckerDelta[m, -1])/2, {m, -nn, nn,
1}];
var = Table[f[i], {i, -nn, nn, 1}];
soli = NDSolve[{eq, ic}, var, {t, 0, tm}]; // AbsoluteTiming


Here is a comparison with the calculation using code @user55777 (left), Fourier (top right) and both codes at time t = 30 (bottom left). The bottom right shows the calculation without integrals (green curve) and with integrals (red curve) by the Fourier method.

Plot3D[Evaluate[Re[uf[x, t] /. soli]], {x, -Pi, Pi}, {t, 0, tm},
Mesh -> None, ColorFunction -> "Rainbow"]


• @user55777 (1) uf[x,t] used in Plot3D[Evaluate[Re[uf[x, t] /. soli]], {x, -Pi, Pi}, {t, 0, tm}, Mesh -> None, ColorFunction -> "Rainbow"]; (2) just call FourierCoefficient[Cot[(x - xp)/2], xp, n] Jul 25 '19 at 8:10
• (1) In that case, @Alex Trounev why you used (1 - Sign[n]) instead of (1 + Sign[n]) in the code? (2) if the computation domain in your code is [-\pi,\pi] how did you plot for [-30,30]? Jul 25 '19 at 9:07
• @user55777 (1) This is due to the definition of the function uf[x_, t_] := Sum[f[n][t] Exp[I n x], {n, -nn, nn}], thus, it is necessary to make a replacement n->-n to satisfy the definition FourierCoefficient[Cot[(x - xp)/2], xp, n]. (2) just use Plot[{Re[uf[Pi x/L, tm] /. soli], u[x, tmax] /. sol}, {x, -L, L}, PlotLegends -> {"Fourier", " @user55777"}]. Jul 25 '19 at 9:35
• @user55777 (1) $u^2u_x$=I Sum[ If[Abs[m - k1 - k2] <= nn, f[k1][t] f[k2][t] (m - k1 - k2) f[m - k1 - k2][t],0], {k1, -nn, nn, 1}, {k2, -nn, nn, 1}]; (2) $u^5(u_x)^2$=- Sum[ If[Abs[m - k1 - k2-k3-k4-k5-k6] <= nn, f[k1][t] f[k2][t] f[k3][t] f[k4][t] f[k5][t] f[k6][t] k6 (m - k1 - k2-k3-k4-k5-k6) f[m - k1 - k2-k3-k4-k5-k6][t],0], {k1, -nn, nn, 1}, {k2, -nn, nn, 1}, {k3, -nn, nn, 1}, {k4, -nn, nn, 1}, {k5, -nn, nn, 1}, {k6, -nn, nn, 1}] Jul 26 '19 at 4:31
• @user55777 Yes, it is right. Jul 26 '19 at 5:12