Troubles in FindFit and NMinimize with NIntegrate

Question

In my code, I have the main function with NIntegrate. When I used Findfit to determine a, b, c and d, there are some error messages show up like

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.

I understand that there may be some singularities in calculations because the point will change at very large R (about = 10^50). I can't remove it by just setting the Exclusions. It's also difficult to converge...

The following is my code, Do you have any idea for it to run?

dat1 = {{48.1059, 0.01599}, {50.5671, 0.00873}, {52.5155,0.00152}, {54.7716,-0.00911}, {56.7201, -0.01437}, {57.7456,-0.01515}, {59.1813, -0.01475}, {60.82214, -0.01125}, {62.4629,-0.00587}, {63.5909, -0.00011}, {63.898645,0.00147}, {65.5394,0.00852}, {68.3083, 0.01671}, {69.128704, 0.01838}, {69.7440,0.01859}, {71.0771, 0.01817}, {71.589908, 0.01709}, {73.0256,0.01355}, {75.7944, 0.00278}, {76.6148, -6*10^-5}};
Q = 5.8*10^-5;
B = 20;
rc = 0.05;
rrw = 0.05;
w = (2*Pi)/24;

FL[r_, z_, t_, Kd_, f_, g_, s_, R_?NumericQ, m_?NumericQ] :=
Block[{x, x1, solx},
solx = FindRoot[Tan[x1] == (s*(-R - Kd*x1^2))/(Kd*x1), {x1, (2 m + 1) Pi/2 + 10^-8}];
x = x1 /. solx;
(((Kd^2 x^2 + (-R - Kd*x^2)^2 s^2) BesselK[0, r Sqrt[R] I] Sin[
   x] (g Cos[t*w] + (-R - Kd*x^2) Sin[t*w]))/(x ((-R - Kd*x^2)^2 +
     g) (Kd^2 x^2 + Kd (-R - Kd*x^2) s + (-R - Kd*x^2)^2 s^2) ((-R - 
       Kd*x^2) f BesselK[0, Sqrt[R] I] + Sqrt[R] I BesselK[1, Sqrt[R] I])))*Cos[x z]]

F[r_, z_, t_, a_?NumericQ, b_?NumericQ, c_?NumericQ, d_?NumericQ] := 
Im[(Q/(2*Pi*B*a))*(2 I)/Pi*ParallelSum[NIntegrate[
  FL[r, z, t, (b/a)/(B/rrw)^2, rc^2/(2*rrw^2*c*B), (c*rrw*rrw)/a*w, d/(c*B), R, m], {R, 0,Infinity}, 
PrecisionGoal -> 6 , MaxRecursion -> 20,Method -> {GlobalAdaptive, MaxErrorIncreases -> 10000}], {m, 0,16, 1}]];

Timing[FindFit[dat1, F[10^-3, 0.5, (t + 3.5), a, b, c, d],
 {{a, 10^-4}, {b, 10^-5}, {c, 10^-5}, {d, 0.1}}, t, Method -> "LevenbergMarquardt", Gradient -> "FiniteDifference"]]

NMinimize[{(Sum[(dat1[[i, 2]] - F[10^-3, 0.5, (dat1[[i, 1]] + 3.5), a, b, c, d])^2,
 {i, 1, Length[dat1]}]/(Length[dat1] - 4))^(1/2), {10^-6 < a < 10^-2, 0 < b < a, 10^-6 < c < 10^-4, 10^-2 < d < 10^-1}}, {a, b, c, d}, Method -> {"SimulatedAnnealing"}]

Thank you!

Answer 1

While I assume there must be some theoretical reason for the model you've coded, what if you only need to make predictions rather than interpret the coefficient estimates? If so, the following is much quicker. A simple sine wave seems to provide an adequate fit using the same number of parameters:

dat1 = {{48.1059, 0.01599}, {50.5671, 0.00873}, {52.5155, 0.00152},
  {54.7716, -0.00911}, {56.7201, -0.01437}, {57.7456, -0.01515},
  {59.1813, -0.01475}, {60.82214, -0.01125}, {62.4629, -0.00587}, 
  {63.5909, -0.00011}, {63.898645, 0.00147}, {65.5394, 0.00852}, 
  {68.3083, 0.01671}, {69.128704, 0.01838}, {69.7440, 0.01859}, 
  {71.0771, 0.01817}, {71.589908, 0.01709}, {73.0256, 0.01355}, 
  {75.7944, 0.00278}, {76.6148, -6*10^-5}};

nlm = NonlinearModelFit[dat1, a + b Sin[2 π x/c + d], 
  {{a, 0.002}, {b, 0.016}, {c, 24}, {d, 2}}, x];
solution = nlm["BestFitParameters"]
(* {a -> 0.0018061570789842571, b -> 0.01691787655549368, 
    c -> 23.841890576175945, d -> 1.954400100333363} *)
Show[ListPlot[dat1], Plot[a + b Sin[2 π x/c + d] /. solution, {x, 48, 76.6}]]