0

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] == 0, Phi1[0] == 0, Phi2[0] == 0, 
  Phi3[0] == 8*Pi*HeteroNuclDen[Tiso], 
  xi[0] == 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.

Luigi
  • 1,301
  • 8
  • 14
  • 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
  • it is intended to be like this (Phi3' and cascading structure). It describes the growth process of randomly formed nuclei (during crystallization). I am not sure how to get another model. Is there a way to get around the underflow? – Luigi May 24 '19 at 11:24
  • You could look for an adequate scaling of time t and functions Phi#. – Ulrich Neumann May 24 '19 at 12:17

1 Answers1

1

Try ParametricNDSolveValue 

Xi = ParametricNDSolveValue[{eq1, eq2, eq3, eq4, eq5, Phi0[0] == 0, Phi1[0] == 0, Phi2[0] == 0, Phi3[0] == 8*Pi*HeteroNuclDen[Tiso],xi[0] == 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 !

Ulrich Neumann
  • 53,729
  • 2
  • 23
  • 55
  • 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]==0,Phi1[0]==0,Phi2[0]==0,Phi3[0]==8 E^(-constN (19.333 +Times[<<2>>])) Nref \[Pi],xi[0]==0},{},{},{},{}}) is less than WorkingPrecision (25.). >>` – Luigi May 24 '19 at 11:46