What is the reason for failing to solve the following equation?

Question

I have a simple differential equation $$ \frac{df(T)}{dT} = g(T)(f(T)-f_{\text{eq}}(T)), $$ where $g,f_{\text{eq}}$ are some smooth functions, and $T$ is an analog of time, but the evolution corresponds to its decreasing. The initial condition is fixed at some large $T$.

As far as $g(T)$ is large, $f(T)$ tends to $f_{\text{eq}}(T)$.

The data and interpolations for $g, f_{\text{eq}}$ is given below:

datatable = {{0.001, -18.23707929395402, -3.4611075922132485*^51, 
    0.06611490804695141}, 
     {0.0012589254117941675, -17.84461494544366, \
-1.731220410720289*^51, 0.06940434750265481}, 
     {0.001584893192461114, -17.431914718313582, \
-8.51731766948239*^50, 0.07237676798672106}, 
     {0.001995262314968879, -17.005104472877587, \
-4.3165099986682115*^50, 0.07395569051393182}, 
     {0.0025118864315095794, -16.557706734429065, \
-2.1593973167434227*^50, 0.07467984975288536}, 
     {0.0031622776601683794, -16.092497328595197, \
-1.074810138034652*^50, 0.07614210743368215}, 
     {0.003981071705534973, -15.612016450430032, \
-5.402227503991311*^49, 0.07670151429807713}, 
     {0.005011872336272725, -15.11979376331143, \
-2.7075274581505683*^49, 0.0770387920958333}, 
     {0.00630957344480193, -14.619694946405195, \
-1.3569781967203665*^49, 0.07725366322862573}, 
     {0.007943282347242814, -14.115080888527912, \
-6.801001485068046*^48, 0.07739012922334851}, 
     {0.01, -13.607560271569687, -3.3982958779998454*^48, 
    0.0771018106537571}, 
     {0.012589254117941675, -13.095519359896302, \
-1.710097520600141*^48, 0.07686997193087373}, 
     {0.01584893192461114, -12.57502890910125, \
-8.690391178966214*^47, 0.0759855227210446}, 
     {0.0199526231496888, -12.044506899327656, \
-4.371902972077502*^47, 0.07282260228552294}, 
     {0.025118864315095808, -11.508867694748597, \
-2.1773382055232237*^47, 0.06994033587672496}, 
     {0.03162277660168379, -10.974125604097805, \
-1.0633537829077974*^47, 0.06503886133118436}, 
     {0.039810717055349734, -10.444363631208544, \
-5.112432708696619*^46, 0.06099002866348603}, 
     {0.05011872336272725, -9.921345251439233, \
-2.501552056154496*^46, 0.05697924054878366}, 
     {0.06309573444801933, -9.405008535108163, \
-1.1914528506215466*^46, 0.05377917472899783}, 
     {0.07943282347242818, -8.893455863058334, \
-5.952272541408941*^45, 0.05042561999840631}, 
     {0.1, -8.385806548209892, -2.914976254908381*^45, 
    0.04632976895359}, 
     {0.12589254117941676, -7.880454160420435, \
-1.4737468497974092*^45, 0.039566117973202475}, 
     {0.15848931924611143, -7.138243474959937, \
-7.174174780329404*^44, 0.02942287460063761}, 
     {0.19952623149688808, -6.626336945176083, \
-2.965865068779192*^44, 0.017144010394089228}, 
     {0.25118864315095824, -6.108043591225297, \
-9.038610774285183*^43, 0.01313376026977404}, 
     {0.31622776601683794, -5.582083828910504, \
-4.548117030242562*^43, 0.012708083189742772}, 
     {0.3981071705534973, -5.052661452093368, \
-2.1806808000591704*^43, 0.012360362229429312}, 
     {0.5011872336272725, -4.524046758491513, -1.112629480001722*^43, 
    0.01199402899466483}, 
     {0.6309573444801936, -3.9994743706180884, \
-5.3097114984886734*^42, 0.011537300885158159}, 
     {0.7943282347242822, -3.479883065855086, \
-2.6955468924524333*^42, 0.011267163642615846}, 
     {1., -2.965140354827245, -1.3165469387026805*^42, 
    0.010824589365892487}, 
     {1.2589254117941675, -2.453828397690545, \
-6.4755820044992164*^41, 0.010584319966540402}, 
     {1.584893192461114, -1.9454167886779412, \
-3.2633903693694033*^41, 0.010215016118577763}, 
     {1.9952623149688808, -1.4391265773321897, \
-1.5547777621410808*^41, 0.010048630681173598}, 
     {2.5118864315095824, -0.9347011259099397, \
-7.989912520361756*^40, 0.009903098920860035}, 
     {3.1622776601683795, -0.4315440818394637, \
-3.8252581674440914*^40, 0.00972893547965847}, 
     {3.981071705534973, 
    0.07068148630045823, -1.9247852842573809*^40, 
    0.00971813690860107}, 
     {5.011872336272725, 0.5721584641714871, -9.626194155274167*^39, 
    0.009606387739884165}, 
     {6.309573444801936, 1.0730475993443855, -4.7648846361866947*^39, 
    0.00964133432146204}, 
     {7.943282347242821, 1.573490988036743, -2.3835091416492644*^39, 
    0.009554013386705201}, 
     {10., 2.073808527167923, -1.1949012290701652*^39, 
    0.009564214098470244}, 
     {12.589254117941687, 2.5742249035666016, -6.027369878562932*^38, 
    0.00955817523693132}, 
     {15.848931924611142, 
    3.0743680427736577, -3.0297492350947725*^38, 
    0.009445529481161313}, 
     {19.952623149688787, 3.574441950151861, -1.4987904366966382*^38, 
    0.009380697335597966}, 
     {25.11886431509582, 4.0744637940244415, -7.670096311434633*^37, 
    0.009133085878059534}, 
     {31.622776601683793, 4.574453067816407, -3.7568252576901522*^37, 
    0.009002785628402312}, 
     {39.81071705534978, 5.074525988991956, -1.898737042484131*^37, 
    0.008740308444495148}, 
     {50.11872336272725, 5.574464995295442, -9.045965296413293*^36, 
    0.008487702523222045}, 
     {63.09573444801943, 6.0745291537589585, -4.569314012352323*^36, 
    0.008328828070587496}, 
     {79.43282347242821, 6.574382102294367, -2.2428280122531448*^36, 
    0.008070902363797764}};

mass = 0.001;
Tini = 79;
Tfinal = 0.001;
rate[T_] = 
  10^Interpolation[datatable[[All, {1, 2}]], InterpolationOrder -> 1][
    T];
jacobian[T_] = 
  Interpolation[datatable[[All, {1, 3}]], InterpolationOrder -> 1][T];
feq[T_] = 
  Interpolation[datatable[[All, {1, 4}]], InterpolationOrder -> 1][T];
fDynamics[f_, T_] = -rate[T]/(1 + 2.2*10^-7 T^6/mass^2)^2 (f - feq[T]);

When trying to solve the equation:

fsol = f /. 
      NDSolve[{D[f[T], T] jacobian[T]^-1 == fDynamics[f[T], T], 
         f[Tini] == 0}, f, {T, Tfinal, Tini}][[1]]

I get the error

NDSolve::ndsz: At T == 79., step size is effectively zero; singularity or stiff system suspected.

Could you please tell me how to avoid it?

Answer 1

The problem is a combination of classic stiffness and weak singularities due to linear interpolation of the data. The stiffness needs a stiff solver plus a higher working precision, as @Ulrich has observed previously and we will explain below. For the weak singularites, NIntegrate has a method "InterpolationPointsSubdivision" to deal with interpolating functions, but NDSolve does not. Such singularities are important because the error estimator assumes the underlying solution is smooth and well-approximated by its Taylor polynomial of the same degree as the order of the method. However, "InterpolationPointsSubdivision" is a simple idea which is easy to implement manually.

First, we increase the precision of the system, which I did initially to avoid having to set and reset the precision; also, it is convenient to have the ODE and initial conditions in separate variables.

{ode, ic} = SetPrecision[
   Rationalize@{
     D[f[T], T] jacobian[T]^-1 == fDynamics[f[T], T], 
     f[Tini] == 0},
   Infinity];

Stiffness and WorkingPrecision

The OP's NDSolve fails out of the gate and cannot even take one step. If we examine the ode at the initial time we can see why. We can also see that MachinePrecision is far from enough precision to handle the problem.

Plot[(f'[T] /. First@Solve[ode, f'[T]]) /. {f[T] -> f} /. 
   Table[{T -> T0}, {T0, 79, 78, -0.25}] // Evaluate
 , {f, 0, 0.01}
 , Frame -> True, FrameLabel -> {f[T], f'[T]}
 , PlotLegends -> Table[T == T0, {T0, 79, 78, -0.25}]]

Since f[T] starts at 0 on the left and, remembering that we are integrating backwards, will rapidly approach but not reach the zero of f'[T] near f[T] == 0.008, the first step must take considerably less time than the following:

dT = 0.008/2.5*^19 (* df/f' *)
(*  3.2*10^-22  *)

dT/T /. T -> 79
(*  4.05063*10^-24  *)

The relative change in T shows that we need well over 24 digits of precision just to be able to take a step at all, but machine precision gives just under 16. It should also be clear that for standard explicit one-step integrators, the step sizes must be small while the derivative stays so large (i.e., the problem will be stiff for such methods).

Interpolation points subdivision for NDSolve

One can use NDSolve`ProcessEquations and NDSolve`Iterate to advance the integration from interpolation point to interpolation point. This will stop the integration on one side of the singular point and restart it on the other. This minimizes the error at the weak singularities. Sometimes, depending on the precision settings and NDSolve method, I would get an error that normally causes NDSolve to abort. The Do loop below would bulldoze right through such errors and usually produce an accurate result. I present settings that do not exhibit such problems.

{state} = NDSolve`ProcessEquations[
   {ode, ic},
   f, {T, Tfinal, Tini}
   , Method -> "LSODA"
   , PrecisionGoal -> 16, AccuracyGoal -> 25, WorkingPrecision -> 50, 
   MaxSteps -> 100000, MaxStepFraction -> 1/3];

stops = SetPrecision[Rest@Reverse@datatable[[All, 1]], 
   state@"WorkingPrecision"];

{newstate} = NDSolve`Reinitialize[state, ic];
Do[
  NDSolve`Iterate[newstate, tt],
  {tt, stops}] // AbsoluteTiming
solIPS = NDSolve`ProcessSolutions[newstate]
f@"Grid" /. solIPS // Dimensions

(*
{0.332851, Null}

{2156, 1}
*)

If we compare with a straight-through integration, we see that the piecewise approach above is faster and the solution has fewer steps:

{newstate} = NDSolve`Reinitialize[state, ic];
NDSolve`Iterate[newstate, Tfinal]; // AbsoluteTiming
solCont = NDSolve`ProcessSolutions[newstate]
f@"Grid" /. solCont // Dimensions

(*
{2.11116, Null}

{14124, 1}
*)

The quality of the solution

We can visualize the steps below (red dots). The initial steps near T == 79 are very small. One can also see that restarting integration begins with small steps.

ListLinePlot[f /. solIPS, PlotRange -> All, Mesh -> All, 
 MeshStyle -> {AbsolutePointSize[1.6], Red}, PlotStyle -> Black, 
 FrameStyle -> Darker@Blue, Frame -> True, GridLines -> {stops, None}]

It is probably better to examine the solution with a log-scale on the horizontal axis (just add ScalingFunctions -> {"Log", Automatic}). We will do that in our comparison with feq[T] below. Note that as we have already observed at machine-precision, the first few steps occur at the same machine-precision value of T. To see the whole plot, we must include and distinguish them, which we accomplish with WorkingPrecision -> 50, Method -> {"BoundaryOffset" -> False}:

Plot[{feq[T], f[T] /. solIPS} // Evaluate
 , Evaluate@Flatten@{T, f@"Domain" /. solIPS}
 , PlotPoints -> 100, ScalingFunctions -> {"Log", Automatic}
 , PlotStyle -> {Directive[AbsoluteThickness[3.], Red, Dashed], 
   Black}, FrameStyle -> Darker@Blue, Frame -> True, 
 GridLines -> {stops, None}
 , WorkingPrecision -> 50, Method -> {"BoundaryOffset" -> False}, 
 PlotRange -> All]

Answer 2

What you call g[T] in your introduction seems to be equal to g[t]==-rate[t]/(1 + 2.2*10^-7 t^6/mass^2)^2 jacobian[t]. This values are of order 10^26!

You need to increase accuracy inside NDSolve dramatically.

Try

F = NDSolveValue[{D[f[T], T] jacobian[T]^-1 == fDynamics[f[T], T], 
f[Tini] == 0}, f, {T, Tfinal, Tini}, WorkingPrecision -> 50] //Quiet

solution

Plot[{F[t], feq[t]}, {t, Tfinal, Tini}, PlotStyle -> {Automatic, Dashed}]

shows F[t]~ feq[t]as expected.