Additional problems with plotting integrated Hankel functions

Question

My goal is to plot a function psi that is setup as

k \[Element] Reals; k > 0;
Ck = HankelH2[1, k]/(HankelH2[1, k] + I HankelH2[0, k])
CGk = Ck (BesselJ[0, k] - I BesselJ[1, k]) + I BesselJ[1, k] //FullSimplify
FGk = Re[CGk]
GGk = Im[CGk]
\[Psi] = 2/\[Pi] NIntegrate[(FGk Cos[k] - GGk Sin[k])/k Sin[k s], {k,0, \[Infinity]}];
\[Psi]plot =  Plot[\[Psi], {s, 0, 20}, PlotRange -> {{0, 20}, {0, 1}},PlotStyle -> {Thickness[.005], Blue}, AxesLabel -> {"s", "\[Psi](s)"}]

The function plots very poorly after about 30 minutes with frightening messages such as:

"NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in k near {k} = {169.7`}. NIntegrate obtained ...error estimates"

The error given in this picture doesn't seem like a lot, but as you can see further down, this function plots poorly compared to a true plot. I would like to increase resolution at the front end of the plot (small values of s).

The function should look more like these approximations:

\[Psi]a1 = 1 - 0.5 Exp[-0.13 s] - 0.5 Exp[-s]
\[Psi]a1plot = Plot[\[Psi]a1, {s, 0, 200}, PlotRange -> {{0, 25}, {0, 1}}, 
  PlotStyle -> {Thickness[.005], Red}, AxesLabel -> {"s", "\[Psi](s)"}]
\[Psi]a2 = (s^2 + s)/(s^2 + 2.82 s + 0.8)
\[Psi]a2plot =  Plot[\[Psi]a2, {s, 0, 200}, PlotRange -> {{0, 25}, {0, 1}}, 
  PlotStyle -> {Thickness[.005], Orange}, AxesLabel -> {"s", "\[Psi](s)"}]
Show[\[Psi]plot, \[Psi]a1plot, \[Psi]a2plot, PlotRange -> {{0, 20}, {0, 1}}, AxesLabel -> {"s", "\[Psi](s)"}]

EDIT:

By reducing the interval of integration as discussed in Julien Kluge's post (thanks!), the "discontinuity" at s = 2 seems to smooth out into an oscillation. I still don't understand why the solution diverges from the correct values (red and orange plots) for values of s < 15.

FINAL EDIT

I believe the function is plotting correctly because I found an error in the original formula. Once corrected, the final solution agrees with the approximate solutions. Changing the integration limits to 100 helped immensely to speed up the calculations.

Answer 1

Okay, this gonna be a long, empirical post so everyone bring your seats in an upright position and fasten your seatbelts.

This is not the solution to your problem but gives new insights to the problem i think.

So whatever we try to convince NIntegrate to converge on $[0,\infty)$ does not help. There are some Methods which do a better job than the automatic selected one but nothing is really able to do it perfectly. But the Oscillatory Methods seem to do a good job. Why that? Well we'll see.

So lets look how our integration-function looks like.

functions=Table[(FGk Cos[k]-GGk Sin[k])/k Sin[k*s],{s,0,5,0.5}];
Plot[Evaluate[functions],{k,0,20},PlotRange->All,PlotTheme->"Detailed"]

Okay, we begin to understand why the Oscillatory-Methods worked better. But hey, no wonder when you see the Cos, Sin and hidden-Bessel functions there.

But we see another thing: all functions vanish pretty fast even on the smaller $[0,20]$-interval compared to $[0,\infty)$. So perhaps we dont need to let NIntegrate decide how far he has to go to approximate to $\infty$.

Lets choose an maximum arbitrary value ($100$) for the interval and $s=1$ and ploting the integrated values in respect to the integration-interval.

testS=1;
t=Table[2/π NIntegrate[(FGk Cos[k]-GGk Sin[k])/k Sin[k*testS],{k,0,end}],{end,1,100}];
Block[{$MaxPrecision=∞},
ListLinePlot[t,PlotRange->All]
]

Year. Exactly what we expected. Our Integral converges pretty fast.

So we need to find a error-estimation for our arbitrary choosen $[0,100]$-interval. So we make an overview-calculation of the max error in $\epsilon=\max_{0\leq x\leq 100}f(x)-\min_{0\leq x\leq 100}f(x)$ in respect to our parameter $s$:

deltaMeasurement[s_?NumericQ]:=Module[{valueList},valueList=ParallelTable[2/π NIntegrate[(FGk Cos[k]-GGk Sin[k])/k Sin[k*s],{k,0,end}],{end,90,100,0.5}];
Differences[MinMax[valueList]][[1]]]
l=deltaMeasurement/@Range[1,20]
Block[{$MaxPrecision=\[Infinity]},ListLinePlot[l,PlotRange->All]]

Okay, there seems to be a finite $s$ where we get the biggest error. It seems to lay in the neighbourhood of $s\approx2$. So we look again.

deltaMeasurement[s_?NumericQ]:=Module[{valueList},valueList=ParallelTable[2/π NIntegrate[(FGk Cos[k]-GGk Sin[k])/k Sin[k*s],{k,0,end}],{end,90,100,0.5}];
Differences[MinMax[valueList]][[1]]]
l=deltaMeasurement/@Range[1.8,2.2,0.025];
Block[{$MaxPrecision=\[Infinity]},ListLinePlot[l,PlotRange->All]]

Cool. We didn't calculated it perfectly now but you Graph will lay in the magnitude of $\approx 10^0$ but our error will be at maximum by $\max_s\epsilon\approx 2\cdot10^{-3}$. Lets confirm that:

FindMaximum[deltaMeasurement[y],{y,2}]

{0.00137326,{y->2.04086}}

Okay so we've won. We know: if we set the integrationinterval to just $[0,100]$ we expect a maximum error of $\frac{\pi}{2}\cdot 1.37326\cdot 10^{-3}\approx 0.0022$ thats far below your ploting-precision.

So without any bad feelings, we can write:

ψ[s_?NumericQ]:=2/π NIntegrate[(FGk Cos[k]-GGk Sin[k])/k Sin[k s],{k,0,100}];
ψplot=Plot[ψ[sp],{sp,0,20},PlotRange->{{0,20},{0,1}},PlotStyle->{Thickness[.005],Blue},AxesLabel->{"s","ψ(s)"}]

Interestingly we still get the discontinuity at $s\approx 2$. Is that an error? I don't know. But there is definetly a connection. Perhaps a rash in far right after 100?

You can now unfasten your seatbelts.