# Fitting system of differential rate equations to data

I am trying to fit a set of differential equations to experimental data, without luck so far. These are the data:

data={{0.630957, 0.00015356}, {0.794328, 0.000327116}, {1.,0.000696757}, {1.25893, 0.00148378}, {1.99526, 0.00671661}, {2.51189, 0.0142547}, {3.16228,0.0301237}, {3.98107, 0.0630841}, {6.30957,0.256021},{10., 0.738742}, {12.5893,0.942704}, {15.8489, 0.997739}, {19.9526, 0.999998}, {25.1189,1.}, {31.6228, 1.}};


First some functions are defined:

Gmax = 1*10^-4;
constG = 6.7*10^-3;
TGref = 60;
GRate[T_] := Gmax*Exp[-constG*(T - TGref)^2]

Nref=8.82393;
constN=0.819435;
TNref=8.82393;
HeteroNuclDen[T_] := Nref*Exp[-constN*(T - TNref)]


Then, after reading the post How can I fit a differential equation to experimental data? I came to the following algorithm

Tiso = 20;
T0m=140;

eq1 = Phi3'[t] == 8*Pi*(aNH*(Tiso + 273))*Exp[-(538/(Tiso - T0m))]*Exp[-(constNHom/((Tiso + 273)^4*(Tiso - T0m)^2))];
eq2 = Phi2'[t] == GRate[Tiso]*Phi3[t];
eq3 = Phi1'[t] == GRate[Tiso]*Phi2[t];
eq4 = Phi0'[t] == GRate[Tiso]*Phi1[t];
eq5 = xi'[t] == (1 - xi[t])*Phi0[t];

Clear[aNH, constNHom]
model[aNH_?NumberQ,constNHom_?NumberQ] := (model[aNH, constNHom] = Module[{Phi0, Phi1, Phi2, Phi3, xi, t},
NDSolveValue[{eq1, eq2, eq3, eq4, eq5,
Phi0 == 0, Phi1 == 0, Phi2 == 0,
Phi3 == 8*Pi*HeteroNuclDen[Tiso],
xi == 0},
{Phi0[t], Phi1[t], Phi2[t], Phi3[t], xi[t]},
{t, 1, 100}]])


However, if I run

nlm = NonlinearModelFit[data, model[aNH, constNHom][t], {{aNH, 10^9}, {constNHom, 10^-10}}, t, Method -> {NMinimize, Method -> "DifferentialEvolution"}]


I get only errors.

• Your model returns 5 functions wheras your data expects one! – Ulrich Neumann May 24 '19 at 10:18
• ok. how to get all the functions out? – Luigi May 24 '19 at 10:33
• Probably your data describes some points of a function[t], which should be approximated by a function Phi0[t] or Phi1[t] or Phi2[t]or Phi3[t] or xi[t] ??? – Ulrich Neumann May 24 '19 at 10:38
• my goal is that the routine finds the parameters for which xi[t] fits the data – Luigi May 24 '19 at 10:57
• Checking your odes I have some questions: Is Phi3'[t]==const intended? The system of odes shows a cascading structure with a very small factor 2.209851823231375 10^-9 in each equation, which might cause underflow. Perhaps you could provide one ode in xi[t]? – Ulrich Neumann May 24 '19 at 11:17

Try ParametricNDSolveValue

Xi = ParametricNDSolveValue[{eq1, eq2, eq3, eq4, eq5, Phi0 == 0, Phi1 == 0, Phi2 == 0, Phi3 == 8*Pi*HeteroNuclDen[Tiso],xi == 0}, xi  , {t, 1, 100}, {aNH, constNHom}, WorkingPrecision -> 25]


Fitting data to xi[t]

nlm = NonlinearModelFit[data,Xi[aNH, constNHom][t], {{aNH, 10^9}, {constNHom, 10^-10}}, t,Method -> {NMinimize, Method -> "DifferentialEvolution"}]


But the fit is very poor. Probably your ode's are critical !

• this is what I get then: ParametricNDSolveValue::precw: The precision of the differential equation ({{(Phi3^\[Prime])[t$38080]==634515. aNH$38078 E^(-9.40398*10^-15 constNH$38079),(Phi2^\[Prime])[t$38080]==1.541*10^-9 Phi3[t$38080],(Phi1^\[Prime])[t$38080]==1.541*10^-9 Phi2[t$38080],(Phi0^\[Prime])[t$38080]==1.541*10^-9 Phi1[t$38080],(xi^\[Prime])[t$38080]==Phi0[t$38080] (1-xi[t$38080]),Phi0==0,Phi1==0,Phi2==0,Phi3==8 E^(-constN (19.333 +Times[<<2>>])) Nref \[Pi],xi==0},{},{},{},{}}) is less than WorkingPrecision (25.). >>` – Luigi May 24 '19 at 11:46