# Fitting to the implicit function

My data is :

data = {{290, 0.012263719}, {240, 0.007675481}, {220,
0.008038809}, {200, 0.008608707}, {170, 0.010805872}, {150,
0.008832903}, {130, 0.009263129}, {120, 0.011290667}, {100,
0.014344114}, {40, 0.025720622}, {30, 0.028876792}, {20,
0.035088327}, {15, 0.042578946}, {4.2, 0.033039767}};


and want to fit it using a function of this form:

a x ((4 Γ^2 + (y^2-R))/(y^2-R)^3)


My only variable is x, y is a function of x and the dependency is given implicitly by:

y == \[Omega]
Sqrt[1 + \[beta] \[Omega]/ y (1/(Exp[y/x] - 1) + 1/2 - y/(2 \[Omega]))]


The parameters that I want to fit are {a, Γ, β,\[Omega],R} and here is my code:

fitfunc[b_?NumericQ, c_?NumericQ, \[Beta]_?NumericQ,
R_?NumericQ, \[Omega]_?NumericQ, \[CapitalGamma]_?NumericQ,
x_?NumericQ] :=
x (10 b + (c/100000) (y^2 - R))/(y^2 - R)^3 /.
FindRoot[
y == \[Omega] Sqrt[1 + \[Beta] \[Omega]/y (1/(Exp[y/x] - 1) + 1/2 - y/2)], {y,
1}];
sol = NonlinearModelFit[
data, {fitfunc[b, c, \[Beta], R, \[Omega], \[CapitalGamma],
x], {b > 0, c > 0, \[Beta] > 0,
R > 0, \[Omega] > 0, \[CapitalGamma] > 0}}, {{b, 8}, {c,
1.3}, {\[Beta], 57}, {\[Omega], 23}, {\[CapitalGamma], 4}}, x];
sol["BestFitParameters"]
(*{b\[Rule]12800,c\[Rule]200000000,\[Beta]\[Rule]57,R \[Rule]3,\
\[Omega] \[Rule]23,\[CapitalGamma]\[Rule]4}*)
(*Original parameters*)
Solve[{a 4 \[CapitalGamma]^2 == 10 b, a == c/100000} /.
sol["BestFitParameters"], {a, \[CapitalGamma]}][[2]]
(*{a\[Rule] 2000,\[CapitalGamma]\[Rule]4}*)

Show[ListPlot[data],
Plot[fitfunc[b, c, \[Beta], R, \[Omega], \[CapitalGamma], x] /.
sol["BestFitParameters"], {x, 1, 300}]]


I know that for this set {a=2000, Γ=4, β=57,\[Omega]=23 ,R=3}} could be one answer ( I found them by changing parameters with hand):

data = {{290, 0.012263719}, {240, 0.007675481}, {220,
0.008038809}, {200, 0.008608707}, {170, 0.010805872}, {150,
0.008832903}, {130, 0.009263129}, {120, 0.011290667}, {100,
0.014344114}, {40, 0.025720622}, {30, 0.028876792}, {20,
0.035088327}, {15, 0.042578946}, {4.2, 0.033039767}};
Show[ListPlot[data],
Plot [{ 2 10^3 x  (4 4^2 + (y^2 - R^2 ))/(y^2 - R^2 )^3 /.
R -> 3} /.
FindRoot[
y == Subscript[\[Omega], 0]
Sqrt[1 +
57 Subscript[\[Omega], 0]/
y (1/(Exp[y/x] - 1) + 1/2 -
y/(2 Subscript[\[Omega], 0]))] /. {Subscript[\[Omega],
0] -> 23}, {y, 1}], {x, 1, 300}], PlotRange -> {0, 0.05}]


I don't know why it couldn't fit and find this set at least.

• The code you've written doesn't line up with the problem you've stated. Where did b and c come from? Also a x ((4 Γ^2 + (y^2-R)/(y^2-R)^3) has mismatched parentheses and contains a potentially degenerate expression (y^2-R)/(y^2-R)^3. – eyorble Oct 18 '17 at 2:59
• Well considering the above plot it should have an answer around this set '{a=2000, Γ=4, β=57,[Omega]=23 ,R=3}}' and I guess it shouldn't be degenerate at least around this set. – Kassik Oct 18 '17 at 4:34
• You've left off R in the NonlinearModelFit statement and there's no $\Gamma$ in the definition of fitfunc. – JimB Oct 18 '17 at 4:48
• Clarification: You've left off R in the parameter and starting value list in the NonlinearModelFit statement and there's no $\Gamma$ in the definition of fitfunc. – JimB Oct 18 '17 at 5:26

I get a dissimilar set of values that fit adequately well by developing an objective function and using NMinimize on it:

fitfunc[a_?NumericQ, \[Beta]_?NumericQ, \[CapitalGamma]_?NumericQ, R_?NumericQ, \[Omega]_?NumericQ, x_?NumericQ] := a x ((4 \[CapitalGamma]^2 + (y^2 - R))/(y^2 - R)^3)
/. FindRoot[y == \[Omega] Sqrt[1 + \[Beta] \[Omega]/y (1/(Exp[y/x] - 1)) + 1/2 - y/(2 \[Omega])], {y, 1}];

of = Total[(Norm@fitfunc[a,\[Beta],\[CapitalGamma],R,\[Omega],#]&/@(Transpose[data][[1]])-Transpose[data][[2]])^2];

mres = NMinimize[{of,a>0,\[Beta]>0,\[CapitalGamma]>0,R>0,\[Omega]>0},{a,\[Beta],\[CapitalGamma],R,\[Omega]}];

Show[ListPlot[data],Plot[fitfunc[a,\[Beta],\[CapitalGamma],R,\[Omega],x]/.mres[[2]],{x,1,300}, PlotRange->{0,0.05}]]


NMinimize has trouble with it too, but at least on my machine it has a much more reasonable fit than NonlinearModelFit gets. The Norm in the objective function is to avoid issues caused by complex numbers appearing during NMinimize's operation.

Edit: Since adding the extra division pointed out in the comments, the objective function above may fail to fit the first few points as well. If that's an issue, consider replacing the entire line with something like:

of = Total[(1000/Transpose[data][[1]])(Norm@fitfunc[a,\[Beta],\[CapitalGamma],R,\[Omega],#]&/@(Transpose[data][[1]])-Transpose[data][[2]])^2];


The (1000/Transpose[data][[1]]) term will weight the accuracy of the terms near the origin higher, and so this process is more likely to achieve a good fit there.

• Well there is a mismatch in your code with mine, – Kassik Oct 18 '17 at 5:53
• Two points: there is a mismatch in the equation in your code with mine, 'FindRoot[y == [Omega] Sqrt[1 + [Beta] [Omega] y (1/(Exp[y/x] - 1)) + 1/2 - y/(2 [Omega])]' should be 'FindRoot[y == [Omega] Sqrt[1 + [Beta] [Omega]/ y ((1/(Exp[y/x] - 1) + 1/2 - y/(2 [Omega]))]'. Second, after performing this revision, the plot becomes disconnected. It doesn't say what are the values of these parameters. – Kassik Oct 18 '17 at 6:05
• The former is an easy fix. The parameters it fits are the second argument of NMinimize, so {a,\[Beta],\[CapitalGamma],R,\[Omega]}. The parameters are stored as mres[[2]] above. – eyorble Oct 18 '17 at 6:08
• @Kassik: Regarding the plot becoming disconnected, adding PlotRange->{0,0.05} to the Plot options helped for me. Plot doesn't normally plot points that would be outside of its graphic, but Show doesn't inform Plot of what the regions of the other graphics are, so Plot still does its own thing. – eyorble Oct 18 '17 at 6:20
• 'PlotRange->{0,0.05}' is not fixing it because for the set of data stored in 'mres[[2]]', the plot is not convergent with 'NMinimize' – Kassik Oct 18 '17 at 6:47