# NDSolve behaves worse than ParametricNDSolve in FindFit

I am trying to fit some model parameters to some data.

Ts = 0.001; duty = 0.2;
Ubat = 14; Ud = 0.6; tEnd = 0.005; Td = 8 10^-7;
pars = {R1 -> 11, L1 -> 0.01};


This is the RL cuircuit model driven by a PWM voltage.

u[t_, d_, T_, Ub_, Ud_] := If[Mod[t, T] < d T, Ub, -Ud]
eqn[R1_, L1_] = {i1'[t] == (u[t] - R1 i1[t])/L1};
ic[i01_] = {i1[0] == i01};


Make some noisy data.

times = N[Range[0, tEnd, Td 20]];
noisyData =
Transpose[{times,sol[times] + RandomVariate[NormalDistribution[0, 0.01], Length[times]]}];
Show[Plot[sol[t], {t, 0, tEnd}], lp = ListPlot[noisyData, PlotRange -> All,  PlotStyle -> Red, PlotLegends -> {"data"}]]


Than I compare the two different approches:

With normal NDSolve and "memoization":

modelRLNDS[R1_?NumberQ,L1_?NumberQ,i01_?NumberQ]:=
modelRLNDS[R1,L1,i01]=(i1/.NDSolve[Flatten[{eqn[R1,L1]/.u[t]->u[t,duty,Ts,Ubat,Ud],ic[i01]}],i1,{t,tEnd},AccuracyGoal >10][[1]])
AbsoluteTiming[sol1=modelRLNDS[R1,L1,0.14]/.pars]


with ParametricNDSolve:

modelRLParaNDS=ParametricNDSolveValue[Flatten[{eqn[R1,L1]/.u[t]->u[t,duty,Ts,Ubat,Ud],ic[i01]}],i1,{t,tEnd},{R1,L1,i01},AccuracyGoal->10];
AbsoluteTiming[sol2=modelRLParaNDS[R1,L1,0.14]/.pars]


Both yield the same result in the same amount of time. But the models behave different in FindFit:

AbsoluteTiming[
Block[{count = 0}, {fit =
FindFit[noisyData,
modelRLNDS[R1, L1, i01][x], {{R1, 1}, {L1, 0.01}, {i01, 0}}, x,
EvaluationMonitor :> count++, AccuracyGoal -> 5], count}]]


AbsoluteTiming[
Block[{count = 0}, {fit2 =
FindFit[noisyData,
modelRLParaNDS[R1, L1, i01][x], {{R1, 1}, {L1, 0.01}, {i01, 0}},
x, EvaluationMonitor :> count++, AccuracyGoal -> 5], count}]]


The solution with the parametric model is not only correct but it is also faster (only 6 function evaluations)!

Show[Plot[{modelRLNDS[R1, L1, i01][x] /. fit,
modelRLParaNDS[R1, L1, i01][x] /. fit2,
modelRLParaNDS[1, 0.1, 0][x]}, {x, 0, tEnd},
PlotRange -> {{0, tEnd}, {0, 0.4}},
PlotLegends -> {"fit1", "fit2", "starting values"}], lp]


Did I do something wrong with the normal NDSolve?

• The problem with the fit is, that your model is evaluated again and again for every point [x]. To speed up the fit you could use "memoization": model[R1_?NumericQ, L1_?NumericQ] := model[R1, L1] = i1 /. NDSolve[ Flatten[{eqn[R1, L1] /. u[t] -> u[t, duty, Ts, Ubat, Ud], ic}], i1, {t, tEnd}][[1]]. But your fit depends strongly on the starting values ... which would be another question
– grbl
Feb 7, 2017 at 9:31
• Guys, please don't vote to close this question, the update of the question is no longer simple, see my answer for more details. BTW @ohmsweetohm I think you can clean up your question a bit. Feb 8, 2017 at 3:45

OK, the Update of the question turns out to be interesting.

Did I do something wrong with the normal NDSolve?

No, you've done nothing wrong, but you should notice the "normal NDSolve" approach and the ParametricNDSolve approach is different. modelRLNDS is completely a black box, while modelRLParaNDS isn't. For example, modelRLParaNDS can be used to "symbolicly" calculate the derivative of R1, L1, i01 without difficulty:

D[modelRLParaNDS[R1, L1, i01][x], R1]
% /. {R1 -> 1, L1 -> 0.01, i01 -> 0}
Plot[%, {x, 0, tEnd}]


while D[modelRLNDS[R1, L1, i01][x], R1] cannot:

D[modelRLNDS[R1, L1, i01][x], R1]
% /. {R1 -> 1, L1 -> 0.01, i01 -> 0}
Plot[%, {x, 0, tEnd}]


It's not clearly mentioned, but the document of Gradient suggests that, when calculating gradient, FindFit performs symbolic differentiation if the model can be symbolically differentiated, and use finite difference if the model is a blackbox, so it's not surprising that the "normal NDSolve" model and the ParametricNDSolve model gives different fitting result. Still, fit can be improved, by setting "DifferenceOrder" -> 2:

AbsoluteTiming[
Block[{count = 0}, {fit =
FindFit[noisyData, modelRLNDS[R1, L1, i01][x], {{R1, 1}, {L1, 0.01}, {i01, 0}}, x,
EvaluationMonitor :> count++,
Gradient -> {"FiniteDifference", "DifferenceOrder" -> 2}], count}]]
(* {0.222677, {{R1 -> 11.0405, L1 -> 0.0100116, i01 -> 0.141926}, 99}} *)

• Thank you for your answer! Would "memoization" of the ParmetricNDSolve model help to improve the performance? Feb 8, 2017 at 18:17
• @OhmSweetOhm I don't think so, because as far as I can tell ParametricFunction already owns memorization ability. (Just try e.g. modelRLParaNDS[1, 0.02, 0]//AbsoluteTiming twice and check the timing. ) Feb 9, 2017 at 2:28