1
$\begingroup$

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!

$\endgroup$
2
  • 2
    $\begingroup$ slwcon is not an error message. It's a warning. Do you get other messages? -- The warning is probably because the integrand is oscillatory and the ?NumericQ in FL[] prevents symbolic analysis. Given that it's necessary b/c of FindRoot, I don't know what to suggest at the moment. You might be able to use NIntegrate`LevinIntegrandReduce (look up in the NIntegrate Rules tutorial in thedocs) on the formula at the end of FL. Since x depends on R and you have Sin[x] etc. in the integrand, it may not be helpful. $\endgroup$
    – Michael E2
    Aug 7, 2017 at 11:38
  • 1
    $\begingroup$ Looking at the predictions using the starting values suggests that whatever parameters influence the amplitude and phase need changing: predicted = Table[{dat1[[i, 1]], F[10^-3, 0.5, (dat1[[i, 1]] + 3.5), a, b, c, d] /. {a -> 10^-4, b -> 10^-5, c -> 10^-5, d -> 0.1}}, {i, Length[dat1]}]; ListPlot[{dat1, predicted}]. $\endgroup$
    – JimB
    Aug 7, 2017 at 13:39

1 Answer 1

1
$\begingroup$

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

Sine wave fit

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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