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
    $\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 '17 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 '17 at 13:39
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, privacy policy and cookie policy

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