I'm trying to fit data($I$ vs $d$) to expressions like:$$ I(d) = \frac{0.013 d^2 \sin^2(\sqrt{1+\frac{0.013 d^2}{\Omega^2}} \Omega)}{1+\frac{0.013 d^2}{\Omega^2} \Omega^2}, $$ where, $\Omega$ is related to $d$ and $w$ through the following transcendental equation:$$ w ~ d \cos\Omega -\Omega=0. $$ Ultimately, my fit parameter is $w$. I can't simply substitute $\Omega$ in $I(d)$ in terms of $w$ and $d$. I guess I could use FindRoot on the first equation and find an $\Omega$ value for each (d, I) data point and fit the resulting {$d$, $\Omega$} data to the second equation to extract $w$. But, this doesn't seem to be a very efficient way. Is it possible to use NonlinearModelFit in this case? How do I accomplish this? Here's some sample data{d,I} for the fit:

idata = {{5.2, 0.1148}, {6.6, 0.228}, {7.98, 0.183}, {9.44, 
   0.379}, {10.38, 0.373}, {12.17, 0.395}, {13.28, 0.337}, {15.16, 
   0.376}, {16.76, 0.282}}

and the expression for $I(d)$:

I[d_] := (0.013 d^2/Sqrt[(1 + 0.013 d^2)/(\[CapitalOmega]^2)] Sin[
     Sqrt[(1 + 0.013 d^2)/(\[CapitalOmega]^2)] \[CapitalOmega]])^2
  • 2
    $\begingroup$ Leaving a comment here since i tried to work on this. So my approach was to solve the omega from a given {d,w} pair of a fit and return the first function value with it. While this should work, I had alot of problems concerning the fact that your transcendental equation allows for multiple solutions and its not clear how to select the right one when presented in this way. One approach I took to solve that was to first estimate guesses from the first equation and then using the nearest solution to that for the respective d but it turns out, this doesn't fit the data really well, unfortunately. $\endgroup$ Aug 18, 2021 at 16:45
  • $\begingroup$ It's best not to use I for a function name; Mathematica uses it for $i = \sqrt{-1}$. $\endgroup$ Aug 18, 2021 at 17:59
  • $\begingroup$ "But, this doesn't seem to be a very efficient way." If you only have 9 data points, "efficiency" can't possibly be a problem. If you had 1,000 data points, then maybe. Your direct original approach should work (other than the issues raised by @JulienKluge). $\endgroup$
    – JimB
    Aug 18, 2021 at 22:45

2 Answers 2


ParametricNDSolve method

You can parlay ParametricNDSolveValue into a curve-fitting function. Saying that $I(d) = f(d)$ for all $d$ is equivalent to saying that $I'(d) = f'(d)$ for all values of $d$ and $I(d_0) = f(d_0)$ for one value of $d_0$. This makes the equation into an ODE. The constraint $w d \cos \Omega - \Omega$ can be interpreted by NDSolve as an algebraic constraint on the variable $\Omega(d)$ that appears in the differential equation, and $w$ is then a parameter of the ODE system.

f[d_] = (0.013 d^2/Sqrt[(1 + 0.013 d^2)/(\[CapitalOmega][d]^2)] Sin[
      Sqrt[(1 + 0.013 d^2)/(\[CapitalOmega][d]^2)] \[CapitalOmega][d]])^2;
intfn = ParametricNDSolveValue[
          {D[int[d], d] == D[f[d], d], int[0] == 0, 
            w d Cos[\[CapitalOmega][d]] - \[CapitalOmega][d] == 0}, 
          int, {d, -20, 40}, {w}]

enter image description here

However, applying NonlinearModelFit to this function causes my kernel to crash. Allowing for an offset in d gives some better (though still not all that good) results:

fit = NonlinearModelFit[idata, intfn[w][d - d0], {{w, 0.05}, {d0, 6}},
 Plot[fit[d], {d, 5, 18}, PlotStyle -> ColorData[97][2]], 
 PlotRange -> All]
Estimate Standard Error t-Statistic P-Value
w 0.028779 0.00113879 25.2714 3.87824e-8
d0 -5.98175 0.374303 -15.981 9.12007e-7

enter image description here

InverseFunction method

$\Omega(d,w)$ can be defined to be the inversion of the function $w = \Omega/d \cos \Omega$. InverseFunction provides a "formal" inverse that Mathematica can evaluate quickly and which can be used in NonlinearModelFit.

Om[d_] = InverseFunction[#/(d Cos[#]) &]
f[d_, w_] := (0.013 d^2/Sqrt[(1 + 0.013 d^2)/(Om[d][w]^2)] Sin[
 Sqrt[(1 + 0.013 d^2)/(Om[d][w]^2)] Om[d][w]])^2
fit = NonlinearModelFit[idata, f[d - d0, w], {{w, 0.05}, {d0, -5}}, d]

This yields the same results for the estimates of $w$ and $d_0$, though oddly it doesn't return the same values for the other statistics and I'm not sure which version to trust.

  • $\begingroup$ Thanks, this is exactly what I was looking for! If I understand this correctly, ` Om[d]` effectively is an operator which acts on each trial ` w` as ` f[d,w]` gets called inside the fit, right? $\endgroup$
    – Devd
    Aug 19, 2021 at 13:40
  • $\begingroup$ @Devd: Yes, that's right. The fancy computation term for it is a pure function or an anonymous function. So Om takes a parameter d and yields a pure function Om[d] which can then be applied to w. $\endgroup$ Aug 19, 2021 at 13:42

One can obtain an estimate for $w$ by minimizing a sum of squares.

idata = {{5.2, 0.1148}, {6.6, 0.228}, {7.98, 0.183}, {9.44, 
    0.379}, {10.38, 0.373},
   {12.17, 0.395}, {13.28, 0.337}, {15.16, 0.376}, {16.76, 0.282}};

(* Prediction for a specif value of d and w *)
f[d_, w_] := Module[{o},
  o = o /. FindRoot[d w Cos[o] - o == 0, {o, 0.5}];
  0.013 d^2 Sin[o Sqrt[1 + 0.013 d^2/o^2]]^2/((1 + 0.013 d^2/o^2)/o^2)]

(* Sum of squares of observed minus predicted *)
sumOfSquares[w_] := Module[{d, s, o, predicted},
  d = idata[[All, 1]];
  o = s /. FindRoot[{# w  Cos[s] - s == 0}, {s, 0.5}] & /@ d;
  predicted = 0.013 d^2 Sin[o Sqrt[1 + 0.013 d^2/o^2]]^2/((1 + 0.013 d^2/o^2)/o^2);
  Total[(idata[[All, 2]] - predicted)^2]]

Show[Plot[f[d, wEstimate], {d, Min[idata[[All, 1]]], Max[idata[[All, 1]]]}, 
  PlotStyle -> Red],

Sum of squares vs w

We see that the value of $w$ that minimizes the sum of squares is around 0.1. That can be made more precise:

wEstimate = w /. FindMinimum[sumOfSquares[w], {{w, 0.1}}][[2]]
(* 0.0992266 *)

A plot of the data and fit suggests the model is probably not very good:

Show[Plot[f[d, wEstimate], {d, Min[idata[[All, 1]]], Max[idata[[All, 1]]]},  PlotStyle -> Red],

Data and fit

  • $\begingroup$ Thanks, this is very helpful. The expression I used is a simplified model with some other parameters being substituted with typical values. I will try the full fit with all the free parameters. $\endgroup$
    – Devd
    Aug 19, 2021 at 13:49
  • $\begingroup$ If all you have is 9 data points, then adding in additional parameters might get you a good reproduction of the data points but the estimated coefficients might not be very good. In other words, you might end up with an overfitting issue. $\endgroup$
    – JimB
    Aug 19, 2021 at 18:13

Your Answer

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

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