# Function with a sharp resonance which Mathematica fails to integrate. Why?

I'm a new user of Mathematica and I'm trying to use it to calculate the collisional cross-section as a function of energy for a given potential and decay rate. I know that the resulting function should be continuous, but when Mathematica does the integration and plots the function, it gives a discontinuous function. Please help! Here is my code:

rc[k_?NumericQ, b_?NumericQ] := Level[NSolve[3/x^12 - 10/x^6 + (5 k/2) b^2/x^2 == k, x,
Reals], {-1}][[-1]]

i2[k_?NumericQ, b_?NumericQ] := NIntegrate[(1/x^12)/Sqrt[N[k] - (3/(N[x])^12 -
10/(N[x])^6 + (5 N[k]/2) (N[b])^2/(N[x])^2)], {x, rc[N[k], N[b]] + 10^-12,Infinity}]

p2[k_?NumericQ, b_?NumericQ] := 1 - E^(-2 i2[N[k], N[b]])

pb2[k_?NumericQ, b_?NumericQ] := N[b]*N[p2[N[k], N[b]]]

s2[k_?NumericQ] := (N[2 Pi])*NIntegrate[pb2[N[k], N[b ]], {b, 0, Infinity},
MinRecursion -> 3, MaxRecursion -> 20]


The typical probability function for a given energy looks like this:

Plot[p2[5, b], {b, 0, 3}]


Notice that there is a spike somewhere near b=1, where the probability jumps up to 1. However, we know that this function is continuous because the probability for the decay of an orbiting resonance is 1, and the probability will continuously approach 1 for energies both slightly above and slightly below the resonance. Mathematically, the exponent -2i2[k,b] just goes to -Infinity and so 1-E^(-2i2[k,b]) goes to 1.

Since all of these functions are continuous in this way, the function pb2[k,b] should also be continuous, having a value of 0 for b=0 (since b approaches 0 and p2[k,b] approaches something between 0 and 1), so we should be able to integrate pb2[k,b] with respect to b numerically without any issues, thus obtaining s2[k] but when I tried to compute s2[k] for various k's and plot the results, this is what I got:

Plot[s2[k], {k, 1, 10}]


Now everything we know about the mathematics and the physics of the situation suggests that the cross-section should be a continuous function of the energy, so what gives? Why did Mathematica integrate a function which is clearly continuous in both of its variables and give me a /dis/continuous function?

I think it must have something to do with that last numerical integration, because I was given the following error message:

NIntegrate::slwcon: Numerical integration converging too slowly; suspect
one of the following: singularity, value of the integration is 0, highly
oscillatory integrand, or WorkingPrecision too small. >>


Like I said, I'm a beginner and don't know how to deal with this stuff. What should I do to get the nice continuous s2[k] that I want?

-
Possibly it is just a loss-of-precision issue, as suggested by the error message. The numerics in your code should be improved to permit higher precision (numeric values should not be wrapped in N, for example, as this resets them to machine precision). – Oleksandr R. Apr 24 '14 at 17:28
I tried taking out all the N's, but it didn't help. The results were essentially the same. Those two regions between 3.5 and 4.5 were still discontinuously separated from the rest of the graph. There's also another region like this around .01. How could I change the settings of NIntegrate to get a smooth function out of this? – Woody Apr 25 '14 at 10:52
@Oleksandr R. Increasing MinRecursion to 5 and MaxRecursion to 50 also didn't help. You can tell from looking at this plot of pb2[4,b]: plot[pb2[4,b],{b,0,3}] ![plot of pb2[4,b]](C:\Users\gsasha\Downloads\Mathematica%2c probability for k%3d4.jpg). Clearly, you can see on this plot that the area under the curve is less than 1/2, probably closer to 1/4, so the value of s2[4] should be less than Pi and closer to Pi/2, but yet when Mathematica computes s2[4] It returns 3.49941, which is greater than Pi, so something is still wrong with the numerical integration... – Woody Apr 25 '14 at 11:13
I also tried taking out the MinRecursion->5, MaxRecursion->50, and it gave me just a smooth decreasing curve without the jumps in the middle. Why should making the integration more accurate make the graph look less continuous? – Woody Apr 27 '14 at 7:23
Your question is a good one, and it is not surprising that as a beginner you have run into trouble here. There are precision issues, but this is not the only problem. I've been looking into it, but with limited success so far due to the extremely long time it takes to produce the plots, especially when working in arbitrary precision. I suspect you will get the right answer (or, at least, one that is wrong for the right reasons) by increasing the precision, but so far I haven't had the patience to confirm this. If I don't make progress, I'll just post the corrected code for you to work on. – Oleksandr R. Apr 27 '14 at 16:56