# Problem with fit on a numerical function

My question is related to a previous one: Problem with simple fit, where I had a problem to fit a simple analytic function $\dfrac{a}{\sqrt{b^2-x^2}}$.

The solution proposed in Problem with simple fit is working very well. However when trying to extend this fitting procedure to another (more complex) list, it does not work.

Here is the list of data I want to fit

listFit = {{10.063521, 0.075566752}, {10.063533,
0.075948022}, {10.063545, 0.076335511}, {10.063557,
0.076620465}, {10.063569, 0.077130577}, {10.063581,
0.077369663}, {10.063594, 0.077928066}, {10.063606,
0.078179029}, {10.063618, 0.078777947}, {10.06363,
0.079075133}, {10.063642, 0.079499388}, {10.063654,
0.079948343}, {10.063667, 0.080276773}, {10.063679,
0.080829161}, {10.063691, 0.08120991}, {10.063703,
0.081737999}, {10.063715, 0.082144721}, {10.063727,
0.082663574}, {10.06374, 0.083099446}, {10.063752,
0.08347671}, {10.063764, 0.084042498}, {10.063776,
0.084695357}, {10.063788, 0.085111893}, {10.0638,
0.085652364}, {10.063813, 0.086138049}, {10.063825,
0.086619555}, {10.063837, 0.087291384}, {10.063849,
0.087728427}, {10.063861, 0.088332497}, {10.063873,
0.088948531}, {10.063886, 0.089623189}, {10.063898,
0.090135636}, {10.06391, 0.09063353}, {10.063922,
0.091418646}, {10.063934, 0.091992204}, {10.063946,
0.092718055}, {10.063958, 0.093453208}, {10.063971,
0.094003733}, {10.063983, 0.094806353}, {10.063995,
0.095456824}, {10.064007, 0.096155563}, {10.064019,
0.096890442}, {10.064031, 0.097742368}, {10.064044,
0.098434713}, {10.064056, 0.099370392}, {10.064068,
0.10007115}, {10.06408, 0.10097375}, {10.064092,
0.10181072}, {10.064104, 0.10272999}, {10.064117,
0.10369551}, {10.064129, 0.1047189}, {10.064141,
0.10558052}, {10.064153, 0.10646513}, {10.064165,
0.10756963}, {10.064177, 0.10864458}, {10.06419,
0.10958263}, {10.064202, 0.11077386}, {10.064214,
0.11191695}, {10.064226, 0.11310339}, {10.064238,
0.11443975}, {10.06425, 0.11560699}, {10.064262,
0.11695356}, {10.064275, 0.11830794}, {10.064287,
0.11967723}, {10.064299, 0.12121895}, {10.064311,
0.122602}, {10.064323, 0.12411718}, {10.064335,
0.12567948}, {10.064348, 0.12756902}, {10.06436,
0.12939693}, {10.064372, 0.13123784}, {10.064384,
0.13290587}, {10.064396, 0.13504858}, {10.064408,
0.13713008}, {10.064421, 0.13917298}, {10.064433,
0.14156279}, {10.064445, 0.14412959}, {10.064457,
0.14645984}, {10.064469, 0.14905608}, {10.064481,
0.15193576}, {10.064494, 0.15481572}, {10.064506,
0.15813123}, {10.064518, 0.16150957}, {10.06453,
0.16513978}, {10.064542, 0.1690908}, {10.064554,
0.1728227}, {10.064566, 0.17751269}, {10.064579,
0.18234043}, {10.064591, 0.18726098}, {10.064603,
0.19313381}, {10.064615, 0.19974816}, {10.064627,
0.20622093}, {10.064639, 0.21374233}, {10.064652,
0.22235096}, {10.064664, 0.23229261}, {10.064676,
0.24326983}, {10.064688, 0.25589852}, {10.0647,
0.2708861}, {10.064712, 0.2892596}, {10.064725,
0.31124154}, {10.064737, 0.33955139}};


Based on Problem with simple fit, I use the following function to fit

fitSingularity2[listFit_] :=
Module[{res, startVal1, sqRootF, model, sf, a, b, fit},
startVal1 = listFit[[-1]][[1]];
res = Sqrt[(startVal1^2 - listFit[[6]][[1]]^2)]*listFit[[6]][[2]];
sqRootF[x_, b_] := b^2 - x^2;
model = a/Sqrt[sqRootF[x, b]];
sf = Reduce[
And @@ Thread[sqRootF[#, b] & /@ listFit[[All, 1]] > 0], {a, b},
Reals];
fit = NonlinearModelFit[
listFit[[2 ;; -1]], {model, sf}, {{a, res}, {b, startVal1}}, x,
Method -> {"NMinimize", Method -> "NelderMead"}][
"BestFitParameters"];
a = a /. fit;
b = b /. fit;
{a, b}];


But as you can see by passing this to Mathematica

fit = fitSingularity2[listFit]
Print[Show[ListPlot[listFit],
Plot[fit[[1]]/Sqrt[((fit[[2]])^2 - x^2)], {x, listFit[[1]][[1]],
listFit[[-1]][[1]]}, PlotStyle -> Red]]]


that the fit is not correct. However I am sure that the function has the form of $\dfrac{a}{\sqrt{b^2-x^2}}$ since by passing this to Mathematica

Print[Show[ListPlot[listFit],
Plot[0.012/Sqrt[((10.06479)^2 - x^2)], {x, listFit[[1]][[1]],
fit[[2]]}, PlotStyle -> Red]]]


I superpose correctly the curves.

• You don't need Print to display results in Mathematica – m_goldberg Sep 22 '15 at 15:04
• Yes you are right. It was previously inside a function it is why the Print is here, I forgot to put it out, my bad. – lambertmular Sep 22 '15 at 16:06

Your x coordinates range are too small and the precision suffers.Just rescale:

list = {#1*100, #2} & @@@ listFit;


Then use my solution to your other question:

sqRootF[x_, b_] := b^2 - x^2
f[x_] := a/Sqrt[sqRootF[x, b]]
sf = Reduce[And @@ Thread[sqRootF[#, b] & /@ list[[All, 1]] > 0], {a, b}, Reals];
nlm = NonlinearModelFit[list, {f[x], sf}, {a, b}, x,
Method -> {"NMinimize", Method -> "NelderMead"}] // Quiet;
nlm["BestFitParameters"]
Show[ListPlot[list],
Plot[nlm[x], {x, Min[list[[All, 1]]], Max@list[[All, 1]]}, PlotStyle -> Red]]

(* {a -> 1.2169239244205383, b -> 1006.4800533971684} *)


Your parameters are obtained by dividing those { a, b } by 100. If in doubt:

Solve[a1/Sqrt[b1^2 - (100 x)^2] == a/Sqrt[b^2 - (x)^2] /.  b -> b1/100, {a}]
(*
{{a -> a1/100}}
*)

• Thank you very much. However, there is one thing that is not clear to me: in the case discussed here the scale is too small and I understand it well. However, for the other kind of singularity ($\dfrac{a}{\sqrt{x^2-b^2}}$) with the same very small scale, a FindFit without any constraint is doing the job perfectly (see fitSingularity function in mathematica.stackexchange.com/questions/95156/…). How is it possible?? – lambertmular Sep 22 '15 at 16:15
• @lambertmular. All software packages can have the same issues when the first 5 digits of the predictor variable are identical when fitting such models. Scaling and centering can many times resolve numerical instability issues. This is not a Mathematica issue. – JimB Sep 22 '15 at 18:30
• @lambertmular. While the fit looks good by eye, there is still a great deal of patter in the residuals: resid = nlm["FitResiduals"]; pred = nlm["PredictedResponse"]; ListPlot[Table[{pred[[i]], resid[[i]]}, {i, Length[pred]}], AxesLabel -> {"Predicted", "Residual"}]. In other words, there seems to be some thing more than the proposed functional form. – JimB Sep 22 '15 at 21:13
• @Jim Baldwin. That's ok here because what I want to do is to fit the function and calculate the integral of the function analytically close to the singularity. This fit should give a good approximation of the integral. But maybe now, thanks to you, I understand why Nintegrate returns 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. when I try to integrate the remaining part of the function. It should be because of the residuals. – lambertmular Sep 23 '15 at 8:58

Following Belasarius's procedure to locate the constraints one finds that b must be > 10.0647.

Plot the error as a function of a and b near the solution.

inco[a_, b_] := Total[Map[(#[[2]] - a/Sqrt[b^2 - #[[1]]^2])^2 &, listFit]]

Plot[{inco[0.01, b], inco[0.011, b], inco[0.012, b], inco[0.013, b]}, {b, 10.0647, 10.07}, AxesOrigin -> {10.064, 0}]

Now blow it up

Plot[{inco[0.01, b], inco[0.011, b], inco[0.012, b],
inco[0.013, b]}, {b, 10.0647, 10.07}, AxesOrigin -> {10.0647, 0},
PlotRange -> {{10.0647, 10.065}, {0, 0.2}},
Ticks -> {{10.0647, 10.0648, 10.0649, 10.065}, Automatic}]


The answer appears to be near 0.012 and 10.0648.

Note also that the error (or incoherence) grows rapidly for values less than 10.0648.

Therefore it looks like a good idea to limit the solution space to a number slightly greater than 10.0647 (for numerical stability purposes).

With these values for constraint and starting value a good solution is derived.

FindFit[listFit, {f[x], b > 10.06475}, {{a, 0.013}, {b, 10.065}}, x,
Method -> {"NMinimize", Method -> "NelderMead"}]

{a -> 0.0121688, b -> 10.0648}


Rescaling

As per Belasarius's solution and Jim Baldwin's comment, rescaling helps with this problem.

listFitR = {#1*100, #2} & @@@ listFit;


Get the constraints

sf = Reduce[
And @@ Thread[sqRootF[#, b] & /@ listFitR[[All, 1]] > 0], {a, b},
Reals]

b < -1006.47 || b > 1006.47


Use FindFit to get the results.

FindFit[listFitR, {f[x], sf}, {a, b}, x,
Method -> {"NMinimize", Method -> "NelderMead"}]


gives the same results as NonLinearModelFit.

{a -> 1.21692, b -> 1006.48}


For this expression each parameter has to be divided by 100 to go back to the original scale.

With[
{
xmin = Min[listFit[[All, 1]]],
xmax = Max[listFit[[All, 1]]]
},
Show[
ListPlot[listFit, PlotRange -> {{xmin, xmax}, Automatic}],
Plot[Evaluate[f[x] /. {a -> 0.0121692, b -> 10.0648}],
{x, xmin, xmax}, PlotStyle -> Red]
]
]
`

• Thank you, but I need to have it automatic to treat a lot of curves that possess the same singularity. Your solution seems good but not suitable for an automatic treatment. Thank you anyway – lambertmular Sep 22 '15 at 16:04