# Multiple NIntegral with singularity

I need to compute the following integral FF[0]. When evaluating FF[0] Mathematica simply takes too much time. The problem may arise due to singularity in the limit 1. I tried specifing different options of NIntegrate, but no luck. Any suggestions why the computation of FF[0] stalls ? I am working on the code below:

FF[z_] := NIntegrate[f1[v] NIntegrate[f4[u], {u, 1, v}], {v, 1, z}];


It involves the functions below:

f1[z_] := z^2/((1 + z + z^2) Sqrt[1 - z^3]);

f2[z_] := (-2 - 3 z^2 + 5 z^3 + 9/2 z^2 Sqrt[1 - z^3]
NIntegrate[f1[u], {u, 1, z}])/(-1 + z^3);

f3[z_] := (z^2 (8 - 12 z - 2 z^3 + 6 z^4 + 3 Sqrt[1 - z^3] (-4 + z^3)
NIntegrate[f1[u], {u, 1, z}]))/(-1 + z^3)^2;

f4[z_] := 1/((1 + z + z^2) Sqrt[1 - z^3]) (16 (-1 + z) (1 + z + z^2) Sqrt[1 - z^3] +
6 (-4 - 4 z - 4 z^2 + z^3 + z^4 + z^5) NIntegrate[f1[u], {u, 1, z}] +
Sqrt[1 - z^3] (-4 (-4 + 2 z + 2 z^2 + 3 z^3) +
4 (1 + z + z^2) NIntegrate[f2[v], {v, 1, z}] -
3 (1 + z + z^2) NIntegrate[f3[s], {s, 1, z}]));

-
There are simple mistakes - For example, NIntegrate takes only numerical values as bounds, but you are trying to pass an analytic bound. You should read the documentation centre carefully. –  Sektor Mar 5 '14 at 22:49
In FF, you do not write NIntegrate twice, that is the wrong syntax. Again, look at the documentation. –  Tom Dickens Mar 6 '14 at 4:56
Thank you. To Sektor: But I do not need the analytic result. I just need FF[0]. –  chrestun Mar 6 '14 at 6:35
To Tom Dickens: If ff1 etc. are simple functions, we can get FF[0] easily, for instance ff1[z_]:=z;ff4[z_]:=z; FF[0]= 0.125. Why do you say we can not write NIntegrate twice? –  chrestun Mar 6 '14 at 6:43
The syntax is NIntegrate[ f1[v] f4[u] ,{v,1,z},{u,1,v}] . Read the manual. –  Tom Dickens Mar 6 '14 at 23:47

One problem is using NumericQ - see What are the most common pitfalls awaiting new users?

Another problem is how many times NIntegrate is being called. To do the 2D integral in FF, f4 will be called tens of thousands of times. Computing f4 in turns does three 1D integrations, which will call f1, f2, f3 several hundred times. And there's still one more layer: f2 and f3 integrate f1.

Finally, to estimate the integral: For singularities of this type -- essentially one over a square root -- there are several possible tricks to eliminate or integrate out the singularity: substitution, Taylor series expansion, or integration by parts. I'll use a partial Taylor series expansion.

Another trick is to use NDSolve instead of NIntegrate. My experience shows that usually NIntegrate will be more accurate for less effort, but in this case, the nested integrations can be performed much more rapidly with a single call to NDSolve.

My guess is that this integral is zero or nearly zero. Tickling the settings makes the result dance around zero. See examples below. The problem is too complicated to tempt me to investigate further.

Suppose the integrand is of the form $f(z)/\sqrt{1-z^3}$. Let $p(x)$ be the Taylor polynomial of $f(x)$ centered at $x=1$: $$p(x) = f(1) + f'(1)\,(z-1) + {1 \over 2!}\, f''(1)\,(z-1)^2 + \cdots + {1 \over n!}\, f^{(n)}(1) \,(z-1)^n$$

We can write the function in the form $${f(x) \over \sqrt{1-z^3}} = {1 \over \sqrt{1-z^3}} p(x) - {1 \over \sqrt{1-z^3}} \left(f(x) - p(x) \right)$$ The first term may be integrated exactly. In the second, the remainder $f(x)-p(x)$ vanishes at $x=1$ with order $n+1$ and the factor $\sqrt{1-z}$ in the denominator may be eliminated.

### Code

The function calcFF0[approxorder] splits the integrand into two parts as above, one using the Taylor polynomial of order approxorder for the rational function factor and the other being the remainder. The first is called f1partialsum and its integral intpartialsum. The remainder is called f1tail and its integral is represented by inttail. The symbols f2i, f3i, f4i represent the integrals of f2, f3, f4. The OP's integrals are translated into the differential equations in diffeqs0; these are then processed to replace f1 by the partial sum and tail.

ClearAll[calcFF0, f1, FF];
f1[z_?NumericQ] := z^2/((1 + z + z^2) Sqrt[1 - z^3]);
FF[z_?NumericQ] :=
NIntegrate[f1[v] f4int[v], {v, 1, z}, Exclusions -> {v == 1},
AccuracyGoal -> 20, PrecisionGoal -> 20, WorkingPrecision -> 50];

calcFF0timings = {};
calcFF0[approxorder_Integer] /; approxorder >= 1 :=
Block[{diffeqs, ivals, fli, zf1i, f2i, f3i, f4i, f1partialsum,
f1tail, intpartialsum, inttail f1sqrt, zf1i0, zf1i1, dzf1i0, dzf1i1,
diffeqs0, z},
calcFF0timings = {};

{f1partialsum, f1tail} =
1/( Sqrt[1 - z] Sqrt[1 + z + z^2]) {#, z^2/(1 + z + z^2)  - #} &@
Normal@Series[z^2/(1 + z + z^2), {z, 1, approxorder}] // Simplify;
AppendTo[calcFF0timings,
"intpartialsum" ->
First@AbsoluteTiming[
intpartialsum =
Assuming[0 < z < 1,
Integrate[f1partialsum /. z -> u, {u, 1, z}]];]
];
Off[Limit::ztest, Limit::ztest1];
f1sqrt =
Limit[(D[intpartialsum, z] - f1tail)/D[Sqrt[1 - z^3], z], z -> 1];
zf1i0 =
Limit[(intpartialsum +
NIntegrate[f1tail /. z -> u, {u, 1, 0}, AccuracyGoal -> 60,
PrecisionGoal -> 60, WorkingPrecision -> 100] -
f1sqrt Sqrt[1 - z^3])/(Sqrt[1 - z^3] (1 - z)), z -> 0];
zf1i1 =
Limit[(intpartialsum - f1sqrt Sqrt[1 - z^3])/(Sqrt[
1 - z^3] (1 - z)), z -> 1];
dzf1i0 =
Limit[D[(intpartialsum + inttail[z] - f1sqrt Sqrt[1 - z^3])/(
Sqrt[1 - z^3] (1 - z)), z] /. {inttail[z] -> zf1i0,
inttail'[z] -> f1tail}, z -> 0];
dzf1i1 =
Limit[D[(intpartialsum + inttail[z] - f1sqrt Sqrt[1 - z^3])/(
Sqrt[1 - z^3] (1 - z)), z] /. {inttail[z] -> 0, inttail'[z] -> f1tail},
z -> 1];
On[Limit::ztest, Limit::ztest1];
ivals = {
f4i[1] == 0,
f3i[1] == 0,
f2i[1] == 0,
inttail[1] == 0,
zf1i[1] == zf1i1
};
diffeqs0 = {
f4i'[z] ==
1/((1 + z + z^2) Sqrt[1 - z^3]) (16 (-1 + z) (1 + z + z^2) Sqrt[
1 - z^3] + 6 (-4 - 4 z - 4 z^2 + z^3 + z^4 + z^5) f1i[z] +
Sqrt[1 - z^3] (-4 (-4 + 2 z + 2 z^2 + 3 z^3) +
4 (1 + z + z^2) f2i[z] - 3 (1 + z + z^2) f3i[z])),
f3i'[z] == (z^2 (8 - 12 z - 2 z^3 + 6 z^4 +
3 Sqrt[1 - z^3] (-4 + z^3) f1i[z]))/(-1 + z^3)^2,
f2i'[z] == (-2 - 3 z^2 + 5 z^3 +
9/2 z^2 Sqrt[1 - z^3] f1i[z])/(-1 + z^3),
inttail'[z] == f1tail,
zf1i'[z] ==
PiecewiseExpand@
If[z == 0, dzf1i0,
If[z == 1, dzf1i1,
D[(intpartialsum + inttail[z] - f1sqrt Sqrt[1 - z^3])/(
Sqrt[1 - z^3] (1 - z)), z] /. inttail'[z] -> f1tail]]
};
diffeqs =
Simplify@
Collect[# // Simplify, {zf1i[z]}] & /@ (diffeqs0 /. {f1i[z] ->
Sqrt[1 - z^3] ((1 - z) zf1i[z] + f1sqrt)});
AppendTo[calcFF0timings,
"f4int" ->
First@AbsoluteTiming[{f4int} = f4i /. NDSolve[{diffeqs, ivals},
{f4i},
{z, 0, 1},
MaxSteps -> 25000, AccuracyGoal -> 20, PrecisionGoal -> 20,
WorkingPrecision -> 50, Method -> {"StiffnessSwitching"}];]
];
(AppendTo[calcFF0timings, "FF[0]" -> First@#]; Last@#) &@
AbsoluteTiming@FF[0]
]


### Examples

Calculation is pretty quick up to order 2. Around order 5 the expressions become cumbersome.

Table[calcFF0[i], {i, 10}]
(*
{-1.3244026180261123588915843875040816322466525948705*10^-19,
3.0015066202265498858065794938827534572730220506990*10^-20,
3.0015066202265498858065794938827534572730220506990*10^-20,
1.8703459094818216335723220477167003403290466171147*10^-20,
1.6186535958160423926232560464349953531508440283849*10^-20,
1.1735419707791214544755429836505597298579609121478*10^-20,
4.7867151751302942415829600894502557404630460424603*10^-21,
-3.6094172859114470277845844574784350053347977553360*10^-22,
-3.6094172859114470277845844574784350053347977553360*10^-22,
-1.0996197631289901595882055987785067444207300272497*10^-22}
*)

calcFF0[20]
(*
-3.5412378870111975729395101114276072893513118648442*10^-23
*)

-
Great answer！I am sorry that I have not seen this answer two weeks ago. Thank you so much. –  chrestun Mar 24 '14 at 21:00