1
$\begingroup$

I am trying to use ParametricNDSolveValue with NonlinearModelFit in order to fit a differential model to some data. I am following the documentation here under Applications/Parameter Fitting: http://reference.wolfram.com/language/ref/ParametricNDSolveValue.html

The documentation example works great for me. The problem is when I use my own model and data, I keep getting an error message that seems strange and out of context to me:

testdata={{1, 0.00009548941732563731}, {101, 0.9234156316342383}, {201, 3.169678865642687}, 
  {301, 5.770066790545615}, {401, 8.007752270115398}, {501, 9.63764824871765}, 
  {601, 10.71272247965258}, {701, 11.380612793541024}, {801, 11.761472037816578}, 
  {901, 11.694455123303454}, {1001, 11.627820089304084}, {1101, 11.56156489420518}, 
  {1201, 11.495687246587122}, {1301, 11.430184947943506}, {1401, 11.365055874810842}, 
  {1501, 11.300297901728584}, {1601, 11.235908918183338}, {1701, 11.171886822658932}, 
  {1801, 11.108229526053263}, {1901, 11.044934960543287}, {2001, 10.982001044546381}, 
  {2101, 10.919425710642873}, {2201, 10.857206933289556}, {2301, 10.795342678813919}, 
  {2401, 10.73383091650653}, {2501, 10.672669644642585}, {2601, 10.611856874496073}, 
  {2701, 10.551390611711986}, {2801, 10.49126888326679}, {2901, 10.43148973333629}, 
  {3001, 10.372051202912106}, {3101, 10.312951353605266}, {3201, 10.254188257327602}, 
  {3301, 10.195759986157272}, {3401, 10.137664640143038}, {3501, 10.079900328269481}, 
  {3601, 10.022465159521184}, {3701, 9.965357246277472}, {3801, 9.90857473007318}, 
  {3901, 9.852115766465182}, {4001, 9.7959785110755}, {4101, 9.740161119622067}, 
  {4201, 9.684661766620447}, {4301, 9.629478652137225}, {4401, 9.574609977715335}, 
  {4501, 9.520053943251227}, {4601, 9.465808730076304}, {4701, 9.411872583200307}, 
  {4801, 9.35824376907317}, {4901, 9.304920554144836}, {5001, 9.251901204865241}, 
  {5101, 9.199183987684325}, {5201, 9.146767169052024}, {5301, 9.094649015418282}, 
  {5401, 9.042827793698757}, {5501, 8.991301804445536}, {5601, 8.940069395555765}, 
  {5701, 8.889128918137416}, {5801, 8.838478723298456}, {5901, 8.788117162146852}, 
  {6001, 8.738042585790575}, {6101, 8.688253345337593}, {6201, 8.638747791895874}, 
  {6301, 8.589524279233093}, {6401, 8.54058120545004}, {6501, 8.491917001537683}, 
  {6601, 8.443530099121851}, {6701, 8.395418929828384}, {6801, 8.347581925283112}, 
  {6901, 8.300017517111874}, {7001, 8.252724136940502}, {7101, 8.205700216394835}, 
  {7201, 8.158944194399854}, {7301, 8.112454560198403}, {7401, 8.066229821998448}, 
  {7501, 8.020268488014702}, {7601, 7.974569066461876}, {7701, 7.929130065554683}, 
  {7801, 7.883949993507835}, {7901, 7.839027358536044}, {8001, 7.794360668782631}, 
  {8101, 7.74994844766569}, {8201, 7.705789267384929}, {8301, 7.661881709002078}, 
  {8401, 7.61822435357887}, {8501, 7.574815782177037}, {8601, 7.53165457585831}, 
  {8701, 7.488739315684421}, {8801, 7.446068582717102}, {8901, 7.4036409616776115}, 
  {9001, 7.36145508077641}, {9101, 7.319509578949874}, {9201, 7.277803092077042}, 
  {9301, 7.236334256036951}, {9401, 7.195101706708642}, {9501, 7.1541040839964705}, 
  {9601, 7.113340057147644}, {9701, 7.072808307729679}, {9801, 7.032507517306886}, 
  {9901, 6.992436367443573}}

RL = 50.; 
Lk = 0.875*10^6; 
v0 = 0.06; 
Rsq = 461; 
width = 70; 
Iss = 1.04; 
δ = 1.*^-11; 
Ibn = 12.5; 
Nhs = 1; 
tch = Sqrt[(width*Lk*Iss)/(2*Sqrt[2]*v0*Rsq*Ibn)]; 
Rch = Lk/tch; 
δ = 1.*^-11; 
Clear[τr]
model = 
  ParametricNDSolveValue[
    {Derivative[1][Iread][t] == 
       R[t]*((Ibn - Iread[t])/(tch*Rch)) - Iread[t]/τr, 
     Derivative[1][R][t] == 
       Piecewise[
         {{(Rch/(tch*Ibn)) * 
           (((Ibn - Iread[t])^2 - Iss^2) /
             Sqrt[((Ibn - Iread[t])^2 - Iss^2*0.5 + Abs[(Ibn - Iread[t])^2  - 
                   Iss^2*0.5])/2. + δ]), 
           R[t] > 0}, 
          {0, 
           R[t] <= 0}}], 
      Iread[0] == 0, R[0] == 0.}, Iread, 
    {t, 0, 10000.}, τr]
fittest = FindFit[testdata, model[τr][t], τr, t]

Here is the (first) error it outputs:

ParametricNDSolveValue, nlnum : The function value {5.86554*10^-9, 13.3632, 9.01097*10^-14, -2.9806*10^-18 - 0.0214555 (-3.46095*10^-17 - 3.46095*10^-17 (Abs')[155.709])} is not a list of numbers with dimensions {4} at {t\$1765520, NDSolve`Iread$25\$1[t$1765520], NDSolve`R\$109\$1[t\$1765520], NDSolve`Iread\$25\$2[t\$1765520], NDSolve`R\$73$2[t\$1765520], NDSolve`s\$1765597[t\$1765520], τr\$1765519} = {0.0000307256, 9.01111*10^-14, 0.000410594, 1.38438*10^-18, 0., -1, 1.}. >>

Because the error states "...not a list with dimensions {4}..." I was sure I was making some kind of syntax error, but as far as I can tell, my syntax is the same as in the documentation. Is it just that my model is too complicated and for some reason this is the error message generated?

Thanks for putting up with my messy raw data and model — I tried to create a simpler version for posting here, but wasn't able to reproduce the error unless I used the full model.

$\endgroup$
  • $\begingroup$ Welcome! To make the most of Mma.SE start by taking the tour now. It will help us to help you if you write an excellent question. Edit if improvable, show due diligence, give brief context, include minimal working example of code and data in formatted form. As you receive give back, vote and answer questions, keep the site useful, be kind, correct mistakes and share what you have learned. $\endgroup$ – rhermans Jun 1 '18 at 20:17
  • $\begingroup$ Is Ib supposed to be defined? $\endgroup$ – rhermans Jun 1 '18 at 20:31
1
$\begingroup$

I found two issues. The first and minor one is that Ib is undefined. I am sure that this happened only in the course of copying the code to this site. If I set Ib = 1 in my copy of the code, then at least ParametricNDSolveValue runs without any issues.

When I call FindFit afterwards, it complains and sends me the same error message as yours. Note that Abs' appears in the errror message. The point is that FindFit tries to apply standard gradient-based techniques for finding the parameters with least squares error. But Mathematica does not know what the derivative of Abs should be (basically it was made so because Abs is not holomorphic). Replace Abs by RealAbs and this issue gets settled---only to produce another one (ParametricNDSolveValue::smpf). I guess we have to pass suitable Method settings to NDSolve in order to account for some discontinuities in the ODE but I have no idea how to do that. Maybe somebody else can tell how to fix that...

$\endgroup$
  • $\begingroup$ Thanks very much for your response. You are correct, Ib is defined in my notebook but just got left out with copy/paste. I will fix it. Good find with the Abs - hopefully I can do some more reading on this new error and figure it out. (or someone with more knowledge than me will come along with some ideas :) ) $\endgroup$ – Kaquel Jun 4 '18 at 13:43
  • $\begingroup$ Well, then just add the definition of Ib to the question so that visitors to the site don't get confused. I'll edit my answer, too, for taking that into account. $\endgroup$ – Henrik Schumacher Jun 4 '18 at 13:50
  • $\begingroup$ I'm still running MM 10, so I don't have RealAbs; but here's a thought: Is it possible to define a rule that sends Abs'[x] to Sign[x] throughout? Or is that how RealAbs already behaves? $\endgroup$ – Michael Seifert Jun 4 '18 at 14:06
  • 1
    $\begingroup$ @MichaelSeifert Sure, one can do Abs' = Sign. One can also use Sqrt[x^2] instead of Abs. Rounding off the corner is actually a good idea to help the optimization algorithm. $\endgroup$ – Henrik Schumacher Jun 4 '18 at 14:09
  • 1
    $\begingroup$ @HenrikSchumacher: It appears to work, too. I'm writing up an answer now. $\endgroup$ – Michael Seifert Jun 4 '18 at 14:10
1
$\begingroup$

The problem appears to arise because of the ill-definedness of Abs[x] near $x = 0$, particularly as regards its derivative. We can try to get around this by "rounding off" the cusp of the function $|x|$ at $x = 0$:

eps = 0.01;
fakeAbs[x_] = Sqrt[x^2 + eps^2];

In the limit eps -> 0, the function fakeAbs[x] is equal to Abs[x]. We can therefore try to run this code for various values of eps, and see if the values found by NonLinearModelFit seem to converge to any particular value.

fakemodel = ParametricNDSolveValue[{Derivative[1][Iread][t] == 
    R[t]*((Ibn - Iread[t])/(tch*Rch)) - Iread[t]/\[Tau]r, 
   Derivative[1][R][t] == 
    Piecewise[{{(Rch/(tch*Ib))*(((Ibn - Iread[t])^2 - Iss^2)/
          Sqrt[((Ibn - Iread[t])^2 - Iss^2*0.5 + 
               fakeAbs[(Ibn - Iread[t])^2 - Iss^2*0.5])/
             2. + \[Delta]]), R[t] > 0}, {0, R[t] <= 0}}], 
   Iread[0] == 0, R[0] == 0.}, Iread, {t, 0, 10000.}, \[Tau]r]
fittest = FindFit[testdata, fakemodel[\[Tau]r][t], \[Tau]r, t]

(* {\[Tau]r -> 18.5261} *)

Reducing eps further does not seem to change the value of $\tau_r$. However, it does not return a very good fit, either:

Show[ListPlot[testdata], Plot[{model[\[Tau]r /. fittest][x], fakemodel[\[Tau]r /. fittest][x]}, {x, 0, 10000}]]

InterpolatingFunction::dmval: Input value {0.204286} lies outside the range of data in the interpolating function. Extrapolation will be used. >>

enter image description here

The InterpolatingFunction error comes from invoking the "real" model at this parameter value; the domain of the InterpolatingFunction returned by ParametricNDSolveValue at this value of $\tau_r$ only appears to go out to about $10^{-6}$ or so. It's also possible that I just used the wrong value of Ib; I set it equal to 1, but perhaps this is way off.

$\endgroup$
  • $\begingroup$ Thanks for the response! It seems we are definitely getting closer. I fixed the error in the original post, but to fix your code just change any remaining "Ib" to "Ibn" or set "Ib=12.5". I am actually able to get quite a good looking fit by fixing this mistake and giving the fit a starting point for tau_r: {tau_r, 17500}. (Note this is the actual tau_r I used when making the test data.) However, weirdly Mathematica throws a FindFit::sszero error. edit: nevermind, the fitting algorithm is just giving up and setting tau_r to whatever I choose originally, not actually fitting. Hence error. $\endgroup$ – Kaquel Jun 4 '18 at 15:28
  • $\begingroup$ @Kaquel The sszero as you've probably guessed is about the sum of squares being essentially zero. It doesn't look like your "data" has any measurement error but only a bit of round-off error. I have to believe that any "real" data won't have this minimal amount of error. Might you consider adding a bit more random error in the simulated data? $\endgroup$ – JimB Jun 4 '18 at 20:48
  • $\begingroup$ @JimB thanks for the suggestion. I have since tried adding various measurement errors between 1-10% and find similar results. $\endgroup$ – Kaquel Jun 13 '18 at 16:03
  • $\begingroup$ I have found a few things that help. I found this thread - in particular the answer at the very bottom from Wolfram: mathematica.stackexchange.com/questions/32455/… Comparing their solution to mine, I noticed they gave Mathematica some conditions in the FindFit function. In particular, tau_r>0 seems to help a lot. I am now able to give a starting value of tau_r=15000, and get a fit value of 16999.9. (The "true" value is 17500.) $\endgroup$ – Kaquel Jun 13 '18 at 16:14
1
$\begingroup$

This is not an answer but rather an extended comment.

Using @HenrikSchumacher 's suggestion of Sqrt[x^2] for Abs[x], the following finds the "good" fit:

model = ParametricNDSolveValue[{Derivative[1][Iread][t] == 
     R[t]*((Ibn - Iread[t])/(tch*Rch)) - Iread[t]/τr, 
    Derivative[1][R][t] == 
     Piecewise[{{(Rch/(tch*Ibn))*(((Ibn - Iread[t])^2 - Iss^2)/
           Sqrt[((Ibn - Iread[t])^2 - Iss^2*0.5 + 
                Sqrt[((Ibn - Iread[t])^2 - Iss^2*0.5)^2])/
              2. + δ]), R[t] > 0}, {0, R[t] <= 0}}], 
    Iread[0] == 0, R[0] == 0.}, Iread, {t, 0, 10000.}, τr];

(* Mean square error function *)
mse[τr_] := 
  Mean[(testdata[[All, 2]] - (model[τr][#] & /@ 
        testdata[[All, 1]]))^2];

(* Find value of τr that minimizes mse *)
sol = FindMinimum[mse[τr], {{τr, 17000}}]
(* {4.2337870614665895`*^-16,{τr\[Rule]17499.999888098355`}} *)

(* Plot results *)
Show[ListPlot[testdata, PlotStyle -> {PointSize[0.02], Red}],
 Plot[model[τr /. sol[[2]]][t], {t, 1, 9901}]]

Data and fit

There are warnings: Warnings

That FindMinimum seems to work and FindFit doesn't will hopefully suggest to others as to what the problem is.

(There also seems to be a local minimum at around $\tau_r = 250$ so a good starting value is essential.)

$\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.