0
$\begingroup$

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.

Improve this question
$\endgroup$
7
  • $\begingroup$ Your model returns 5 functions wheras your data expects one! $\endgroup$
    – Ulrich Neumann
    Commented May 24, 2019 at 10:18
  • $\begingroup$ ok. how to get all the functions out? $\endgroup$
    – Luigi
    Commented May 24, 2019 at 10:33
  • $\begingroup$ 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] ??? $\endgroup$
    – Ulrich Neumann
    Commented May 24, 2019 at 10:38
  • $\begingroup$ my goal is that the routine finds the parameters for which xi[t] fits the data $\endgroup$
    – Luigi
    Commented May 24, 2019 at 10:57
  • $\begingroup$ 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]? $\endgroup$
    – Ulrich Neumann
    Commented May 24, 2019 at 11:17

1 Answer 1

Reset to default
1
$\begingroup$

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 !

Improve this answer
$\endgroup$
1
  • $\begingroup$ 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.). >>` $\endgroup$
    – Luigi
    Commented May 24, 2019 at 11:46

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.