# Fitting data points with a function obtained by ParametricNDSolve+Integrate

My question concerns fitting data points to a function that is a solution of a set of parametric differential equations, then integrated over time. More specifically, I have these parametric differential equations:

picosec = 10^-12;
ti = 0;  (*initial time of the simulation*)
tf = 12000*picosec;   (*final time of the simulation*)
Pump[t_, t0_, τ_] = (Sech[(t - t0)/τ])^2;  (*Definition of the excitation pulse*)

eqns = {xw'[t] == ηinj*(Peakpower*Pump[t, 1000 picosec, 10 picosec])/(
2.55*10^-19) - xw[t]/(10*picosec)*(1 - xp[t]/560.63) - xw[t]/
10^-9 (*DIFFERENTIAL EQUATION 1*),
xp'[t] ==
xw[t]/(10*picosec)*(1 - xp[t]/560.63) -
1.37*10^9*Gainfp*(2 yp[t]*(xp[t] - 280.315) + (xp[t])) - xp[t]/
10^-9(*DIFFERENTIAL EQUATION 2*),
yp'[t] ==
1.37*10^9*Gainfp*(2 yp[t]*(xp[t] - 280.315) + (xp[t])) - yp[t]/(
4.07 picosec)   (*DIFFERENTIAL EQUATION 3*),
xw == 0, xp == 0, yp == 0 (*initial conditions*)};

solFIT = ParametricNDSolve[
eqns, {xw, xp, yp}, {t, ti, tf}, {ηinj, Gainfp, Peakpower},MaxSteps -> Infinity, Method -> "StiffnessSwitching"]


Now my experimental data is proportional to the integral of solution yp over time:

PhotonFunc[ηinj_,Gainfp_,Peakpower_,t_]={yp[ηinj,Gainfp,Peakpower][t]}/.solFIT;

FittingFunction[ηinj_?NumericQ,ηcoll_?NumericQ,Gainfp_?NumericQ,Peakpower_?NumericQ]:=ηcoll/(4.07 picosec)*Integrate[PhotonFunc[ηinj,Gainfp,Peakpower,t],{t,ti,tf}];


In other words, my fitting function has 3 free parameters: 2 of them are from the differential equation definitions (ηinj, Gainfp), and one introduced at the integration stage (ηcoll) to take into account the proportionality between data points and theoretical function. Peakpower represents the independent variable of the model, since every data point is taken for a specific value of Peakpower.

My data is:

 FITDATA = {{3.3*^-6, 3.171234*^-7}, {3.4736842105263158*^-6,
3.3900810000000006*^-7}, {3.774736842105263*^-6,
3.678039*^-7}, {4.052631578947368*^-6,
4.066959*^-7}, {4.388421052631579*^-6,
4.1602330000000005*^-7}, {4.677894736842105*^-6,
4.688141*^-7}, {5.025263157894737*^-6,
5.029277000000001*^-7}, {5.268421052631579*^-6,
5.396957000000001*^-7}, {5.789473684210527*^-6,
5.850883*^-7}, {6.0210526315789475*^-6,
6.178703000000001*^-7}, {6.368421052631579*^-6,
6.590776*^-7}, {6.854736842105263*^-6,
6.949222000000001*^-7}, {7.329473684210526*^-6,
7.842048*^-7}, {7.665263157894737*^-6,
8.33995*^-7}, {8.105263157894736*^-6,
8.631146000000001*^-7}, {8.336842105263159*^-6,
9.144675*^-7}, {9.031578947368421*^-6,
1.0331121000000002*^-6}, {9.436842105263158*^-6,
1.1031741*^-6}, {9.9*^-6,
1.1531064*^-6}, {0.000010536842105263158,
1.2948182*^-6}, {0.000010861052631578949,
1.3427578000000001*^-6}, {0.00001144,
1.4727559000000002*^-6}, {0.000012389473684210526,
1.6650955*^-6}, {0.000014473684210526315,
2.1499883*^-6}, {0.000016789473684210526,
2.7969681000000004*^-6}, {0.000020957894736842108,
3.6925992000000005*^-6}, {0.000022347368421052632,
3.9389504*^-6}, {0.000023736842105263157,
4.1618265000000005*^-6}, {0.000025821052631578948,
4.574431400000001*^-6}, {0.000026747368421052634,
4.7048833000000005*^-6}};


To fit the data, I wrote the line:

NonlinearModelFit[FITDATA, {FittingFunction[ηinj,ηcoll, Gainfp, Peakpower], 0 < ηinj < 1, 0 <ηcoll < 1}, {{ηinj,0.6}, {ηcoll,0.000000008}, {Gainfp, 0.3}}, PeakPower, Method -> "NMinimize"]


But it give me several error, e.g.:

NonlinearModelFit::nrnum: The function value 1/2 ((-4.70488*10^-6+FittingFunction[0.6,8.*10^-9,0.3,Peakpower])^2+(-4.57443*10^-6+FittingFunction[0.6,8.*10^-9,0.3,Peakpower])^2+<<26>>+(-3.39008*10^-7+FittingFunction[0.6,8.*10^-9,0.3,Peakpower])^2+(-3.17123*10^-7+FittingFunction[0.6,8.*10^-9,0.3,Peakpower])^2) is not a real number at {[Eta]inj,[Eta]coll,Gainfp} = {0.6,8.*10^-9,0.3}. >>

I tried to change fitting method, i tried to use FindFit, but I did not get any fit. Sometimes the error is related to the evaluation of the integral (in this case Integrate[] expects numbers instead of symbols), and sometimes is related to complex values in the fit functions like the posted one. I tried to implement solution given in other posts with similar issues but they did not work.

Do you know if there is a smarter way to approach such a problem?

• Your FittingFunction is complicated integral, which Mathematica can't do in closed form, so any further computation also can't be done. One way to go is using Manipulate: plot both FittingFunction and ListPlot[FITDATA] on the same plot as functions of Peakpower (horizontal axis) and use ηinj, ηcoll, Gainfp as Controls for Manipulate to match curves as fine as possible by hand. – Alx Jul 19 at 15:37
• Like this: Manipulate[Show[Plot[FittingFunction[ηinj, ηcoll, Gainfp, Peakpower], {Peakpower, 3.3*^-6, 0.000026747368421052634}], ListPlot[FITDATA]]],{{ηinj, 0.6}, 0.1, 1}, {{ηcoll, 0.000000008}, 0, 0.00000001}, {{Gainfp, 0.3}, 0.1,1}] Your starting values are very close to what can be fitted by hand. – Alx Jul 19 at 15:40