0
$\begingroup$

I was fit the data using differential equation(eqns) and pfun

I tried data fitting (I gave the parameter directly..a=22500, b=10^-22, c=2*10^-40 and scaling factor k=1.04 10^-36)

data = {{0, 0}, {3.69 10^12, 50.65}, {5.67 10^12, 
134.875}, {7.05 10^12, 225.275}, {9.03 10^12, 
381.65}, {1.09 10^13, 509.5}, {1.25 10^13, 595.3}, {1.45 10^13, 
815.325}, {1.76 10^13, 1225}, {2.10 10^13, 1624.125}, {2.46 10^13,
 2018.725}, {2.84 10^13, 2488.775}, {3.18 10^13, 
2942.9}, {3.68 10^13, 3630}, {4.39 10^13, 4558.65}, {5.52 10^13, 
5925.925}, {6.45 10^13, 7044.075}, {7.23 10^13, 
7972.2}, {8.18 10^13, 9119.575}, {9.38 10^13, 10545}, {1.06 10^14,
 11749.1}, {1.28 10^14, 13760.475}, {1.42 10^14, 
15055.65}, {1.57 10^14, 16484.475}};

eqns = {n'[t] == 22500 L - 10^-22 n[t]^2 - 2.0 10^-40  n[t]^3 , n[0] == 0};

pfun = ParametricNDSolveValue[eqns, \!\(\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(30\)]\(n[s]^2 \[DifferentialD]s\)\), {t, 0, 30}, {L}]

then plot the data

Show[LogLogPlot[1.04 10^-36  pfun[L], {L, 10^12, 2*10^14},  
AxesLabel -> {l, K \!\(\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(30\)]\(n[l, t]^2 \[DifferentialD]t\)\)}, ImageSize -> Large, 
PlotStyle -> Red], ListLogLogPlot[data, PlotStyle -> Black], 
LogLogPlot[0.45 10^-23*L^2, {L, 10^12, 2*10^14}], LabelStyle -> 
Directive[Black, Bold, Medium]]

enter image description here


after, I want to fit the data (by computer directly)

eqns = {n'[t] == a L - b n[t]^2 - c n[t]^3 , n[0] == 0};
pfun = ParametricNDSolveValue[eqns, \!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(30\)]\(n[
   t]^2 \[DifferentialD]t\)\), {t, 0, 30}, {a, b, c, k}];

fit = FindFit[data, 
k*pfun, {{a, 22500}, {b, 10^-22}, {c, 2.0 10^-40},{k, 1.04 10^-36}}, L, 
Method -> {NMinimize, Method -> "DifferentialEvolution"}] // Quiet

But the Findfit code errors appear..

For get rid of the errors, how to rivise the above code?

If you give me an another fitting code, I'll accept thankfully.

Thank you.

$\endgroup$
2
  • $\begingroup$ What do the errors tell you? $\endgroup$
    – JimB
    Jan 16, 2019 at 4:50
  • $\begingroup$ Your last pfun depends on 5(!) parameters {a,b,c,k,L}. The second argument of FindFitshould be something like k*pfun[a, b, c, k, L] $\endgroup$ Jan 16, 2019 at 7:57

2 Answers 2

3
$\begingroup$

It turns out that the integral evaluation within ParametricNDSolve is slow due to symbolic processing. I removed it from ParametricNDSolve and evaluated the integral with NIntegrate turning SymbolicProcessing->0. It was >100x faster on my machine doing it this way. I also like using non-dimensional scaling to keep values closer to 1.

I use the subscript $d$ to relate parameters with dimension to their dimensionless counterparts as shown below.

$$\begin{gathered} L = \frac{{{L_d}}}{{{L_{\max }}}} \hfill \\ a = \frac{{{a_d}{L_{\max }}{t_{\max }}}}{{{n_0}}} \hfill \\ b = {b_d}{n_0}{t_{\max }} \hfill \\ c = {c_d}n_0^2{t_{\max }} \hfill \\ {n_0} = \frac{1}{{\sqrt {K{t_{\max }}} }} \hfill \\ \end{gathered} $$

Since the data appear fairly evenly spread out on a log-log scale, I will also do the fit on a log-log scale to achieve a better weighting. Here is how I recast the example given in the OP.

data = {{0, 0}, {3.69 10^12, 50.65}, {5.67 10^12, 
    134.875}, {7.05 10^12, 225.275}, {9.03 10^12, 
    381.65}, {1.09 10^13, 509.5}, {1.25 10^13, 595.3}, {1.45 10^13, 
    815.325}, {1.76 10^13, 1225}, {2.10 10^13, 1624.125}, {2.46 10^13,
     2018.725}, {2.84 10^13, 2488.775}, {3.18 10^13, 
    2942.9}, {3.68 10^13, 3630}, {4.39 10^13, 4558.65}, {5.52 10^13, 
    5925.925}, {6.45 10^13, 7044.075}, {7.23 10^13, 
    7972.2}, {8.18 10^13, 9119.575}, {9.38 10^13, 10545}, {1.06 10^14,
     11749.1}, {1.28 10^14, 13760.475}, {1.42 10^14, 
    15055.65}, {1.57 10^14, 16484.475}};
(* Scale data max L value *)
sdata = data;
sdata[[All, 1]] = data[[All, 1]]/Max@data[[All, 1]];
(* Data looks equally spaced on log-log scale *)
(* so let's scale it (Note: remove {0,0}) *)
lldata = N@Log@sdata[[2 ;; -1]];
(* create parameters for non-dim scaling *)
tmax = 30;
Kfac = 1.04 10^-36;
scalefac = 1;
Lmax = Max@data[[All, 1]];
n0 = 1/Sqrt[Kfac tmax];
(* Scaled equations *)
eqns = {n'[t] == 
    tmax/n0 (22500 Lmax L - 10^-22 n0^2 n[t]^2 - 
       2.0 10^-40 n0^3 n[t]^3) , n[0] == 0};
pfun = ParametricNDSolveValue[ eqns, n, {t, 0, 1}, {L}];
f = Function[{l}, 
   Log@NIntegrate[(pfun[Exp[l]][t])^2, {t, 0, 1}, 
     Method -> {Automatic, "SymbolicProcessing" -> 0}], Listable];
Show[Plot[f[l], {l, -4, 0.5},  
  AxesLabel -> {l, K \!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(30\)]\(n[l, 
         t]^2 \[DifferentialD]t\)\)}, ImageSize -> Large, 
  PlotStyle -> Red], ListPlot[lldata, PlotStyle -> Black], 
 Plot[Log[0.45 10^-23*Lmax^2*Exp[l]^2], {l, -4, 0.5}], LabelStyle -> 
  Directive[Black, Bold, Medium]]
(* Show scaled equations *)
eqns // Simplify
(* {(n^\[Prime])[t]\[Equal]591.9441354553654` \
L-0.0005370861555295748` n[t]^2-0.00019230769230769236` n[t]^3,n[0]\
\[Equal]0} *)

OP Example Recast

It looks we have recast the problem to the mimic the original state. Now, we can formulate the problem to fit $a,b,c,sf$ just using FindFit out-of-the-box.

(* Dimensionless Equations *)
eqns = {n'[t] == a L - b n[t]^2 - c n[t]^3 , n[0] == 0};
pfun = ParametricNDSolveValue[ eqns, n, {t, 0, 1}, {a, b, c, L}];
ffull = Function[{a, b, c, sf, l}, 
   sf Log@NIntegrate[(pfun[a, b, c, Exp[l]][t])^2, {t, 0, 1}, 
      Method -> {Automatic, "SymbolicProcessing" -> 0}], Listable];
(* rewrite in more typical model form *)
model[a_, b_, c_, sf_][l_] := ffull[a, b, c, sf, l]
{time, fit} = 
 Quiet@FindFit[lldata, 
    model[a, b, c, sf][
     l], {{a, 500}, {b, 0.005} , {c, 0.0005}, {sf, 1}}, l] // 
  AbsoluteTiming
(*  {3.18,{a->440.,b->0.00719,c->0.000705,sf->1.13} *)
Show[Plot[model[a, b, c, sf][l] /. fit, {l, -4, 0.5},  
  AxesLabel -> {Log@l, Log[sf \!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(1\)]\(n[l, 
          t]^2 \[DifferentialD]t\)\)]}, ImageSize -> Large, 
  PlotStyle -> Red], ListPlot[lldata, PlotStyle -> Black], 
 Plot[Log[0.45 10^-23*Lmax^2*Exp[l]^2], {l, -4, 0.5}], LabelStyle -> 
  Directive[Black, Bold, Medium]]
(* Rescale Parameters Back to Dimensioned Forms *)
nd = {a, b, c, L};
scalefactors = {Lmax tmax/n0, n0 tmax, n0^2 tmax, 1/Lmax};
(nd /. fit) ~Times~(1/scalefactors)
(*  {16700.,1.34*10^(-21),7.33*10^(-40),(1.57*10^(14)) L} *)

Actual Fit

The fit looks pretty good and it took about 3 seconds on my machine.

$\endgroup$
0
$\begingroup$

The modified script below doesn't gives formal errors but keeps a lot of time processing before outputting some messages now linked to numerical problems.

eqns = {n'[t] == a L - b n[t]^2 - c n[t]^3 , n[0] == 0};
pfun = ParametricNDSolveValue[eqns, \!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(30\)]\(n[
   t]^2 \[DifferentialD]t\)\), {t, 0, 30}, {a, b, c, k, L}];

fit = FindFit[data, k*pfun[a, b, c, k, L], 
{{a, 22500}, {b, 10^-22}, {c, 2.0 10^-40}, {k, 1.04 10^-36}}, L, 
Method -> {NMinimize, Method -> "DifferentialEvolution"}] // Quiet
$\endgroup$

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.