Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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
share|improve this question
1  
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

1 Answer 1

up vote 1 down vote accepted

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}]

Plot of the interpolated function

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}} *)
share|improve this answer
    
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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.