# Integrate on a logarithmic grid

I am trying to solve numerically an integral equation, and from a couple of tests I can see that the result has a strong logarithmic divergence. For my specific problem it is very important to have accurate results in the region of very small arguments, but that doesn't work since I am forced to replace the lower limit of integration (zero, that is) with non-zero values (0.0001 or whatever) otherwise the program crashes.

So I though about integrating on a logarithmic grid, but I don't know how that works in Mathematica. (I have searched in the mathematica references, but could not find an answer). Anyone has an idea?

α = 2.85;
g = (Pi/2) α;
Nf := 2;

cs[x_] := 2 ArcCos[x]/Sqrt[1 - x^2];
csh[x_] := 2 ArcCosh[x]/Sqrt[x^2 - 1];
prefB[p_, k_, d_] :=
(p^2 + k^2 (1 - 1/d^2))/Sqrt[p^2 d^4 - 1/4 ((p^2 + k^2) d^2/k - k )^2];
pieceB[k_, d_] :=
If[d B[k]/k^2 < 1, cs[d B[k]/k^2], If[d B[k]/k^2 > 1, csh[d B[k]/k^2], 2]];

B[p_] = p^2 ;
iterstep :=
(values =
Parallelize[
Table[{p, p^2 +  g/(Pi^3  Nf) (NIntegrate[
prefB[p, k, d] ((d^2 B[k]^2/k^4 - 1) (Pi - g pieceB[k, d]) +
B[k ]/k^2 d  g^2 csh[g])/(d^2 B[k]^2/k^4 + g^2 - 1),
{d, 0, 1/(1 + p)}, {k, p d/(d + 1), p d/(1 - d)},
WorkingPrecision -> 16,
PrecisionGoal -> 2,
MaxRecursion -> 100,
AccuracyGoal -> 16,
Method -> {"SymbolicPreprocessing", "OscillatorySelection" -> False}] +
NIntegrate[
prefB[p, k, d] ((d^2 B[k]^2/k^4 - 1) (Pi - g pieceB[k, d]) +
B[k ]/k^2 d  g^2 csh[g])/(d^2 B[k]^2/k^4 + g^2 -1),
{d, 1/(1 - p), Infinity},
{k, p d/(d + 1), p d/(d - 1)},
WorkingPrecision -> 16, PrecisionGoal -> 2,
MaxRecursion -> 100, AccuracyGoal -> 16,
Method -> {"SymbolicPreprocessing","OscillatorySelection" -> False}])},
{p, 0, 0.99999, 1/20}]];
B[p_] = Interpolation[values , p, InterpolationOrder -> 4,   Method -> "Hermite"])

Do[iterstep, {3}] // AbsoluteTiming

-
Can you include some code containing your specific problem? –  geordie Mar 18 '13 at 9:12
If a certain input to NIntegrate causes Mathematica to crash, that should be reported to WRI as a bug. In any case I agree with @geordie: without a concrete example this question is going to be difficult to answer in a way that is likely to help you. –  Oleksandr R. Mar 18 '13 at 10:29
As said it is not just an integral, it is an integral equation. In fact I have uploaded the code the other day mathematica.stackexchange.com/questions/21303/… There I was asking a different question but since nobody answered (the user belisarius said he has an idea, but at the end he has probably forgotten) I have spent some time time thinking of alternative solution and got the idea of working on a logarithmic scale directly withing the integral (but I don't know how to implement this, I don't even know if it is possible...) –  Micha Mar 18 '13 at 10:46
I have just edited the previous question to remove one of the two codes I have uploaded the first time. In the meantime I have managed to implement the convergence with a while loop (so I could solve the problem from the title, so to speak), but the problem of the integration limits remains :( –  Micha Mar 18 '13 at 10:57
What Oleksandr and geordie are saying is "I don't understand your question". Try to give a small example of what you'd like to do (try not to hurl a page of code at us! just construct a minimal example). –  acl Mar 18 '13 at 11:28

Like the others I'm not 100% sure of what the problem is. When I run your code it looks fine. I slightly modified it:

values = ParallelTable[
{p, p^2 + g/(Pi^3 Nf) (NIntegrate[prefB[p, k, d] ((d^2 B[k]^2/k^4 - 1) (Pi - g pieceB[k, d]) +
B[k]/k^2 d g^2 csh[g])/(d^2 B[k]^2/k^4 + g^2 - 1),
{d, 0, 1/(1 + p)}, {k, p d/(d + 1), p d/(1 - d)}] +
NIntegrate[prefB[p, k, d] ((d^2 B[k]^2/k^4 - 1) (Pi - g pieceB[k, d]) +
B[k]/k^2 d g^2 csh[g])/(d^2 B[k]^2/k^4 + g^2 - 1),
{d, 1/(1 - p), Infinity}, {k, p d/(d + 1), p d/(d - 1)}])},
{p,0,.98,.02}];


I get lots of warnings about non-convergence and complex values. I assume that the result is supposed to be real, because the imaginary part is very small. I'm not worried about that too much, since if I chop of the imaginary part and plot the function, it looks pretty regular:

f = Interpolation[Re@values];
Plot[f[x],{x,0,1}]


Now I understand you are interested in the value at p=1, which is where Mathematica has trouble if you simply extend the range of p to 1. But you can get close by replacing {p,0,.98,.02} by this

{p, Range[0, .98, .02]~Join~(1 - 10^-Range[2, 8, .5])}


The plot looks pretty much the same, and by examining the table you can see that it seems to converge to 0.979533.

f[1]
(* => 0.979533 *)
Re@values[[-3;;]]
(* => {{1., 0.979532}, {1., 0.979533}, {1., 0.979533}} *)

-
Volker,thank you very much! It is true, I went around in circles a couple of times but the idea is to make the damn thing work, in a way or another. I didn't know about Join, now I am just running the code (it doesn't give errors yet but I suppose it will produce something). Now, if I may ask you another stupid question, how can I tell the iteration to stop when the convergence is achieved, for example when the absolute value between two consecutive solutions is less then say 10^-3? Again, thank you very much, you are pretty much saving my life :) –  Micha Mar 18 '13 at 16:33
In that case I'd probably use a while loop instead of a table... but then it wouldn't work in parallel. –  Volker Mar 18 '13 at 20:02