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

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?

• Try Method-> "StiffnessSwitching" and NDSolve evaluates without error message. – Ulrich Neumann Dec 29 '20 at 14:41
• @UlrichNeumann : but it gives some shit instead of the solution ($f(T)\to f_{\text{eq}}(T)$). – John Taylor Dec 29 '20 at 14:50
• Sorry I cannot check your ode. Where is the definition of g[t]? – Ulrich Neumann Dec 29 '20 at 15:07

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 NDSolveProcessEquations and NDSolveIterate 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} = NDSolveProcessEquations[
{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} = NDSolveReinitialize[state, ic];
Do[
NDSolveIterate[newstate, tt],
{tt, stops}] // AbsoluteTiming
solIPS = NDSolveProcessSolutions[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} = NDSolveReinitialize[state, ic];
NDSolveIterate[newstate, Tfinal]; // AbsoluteTiming
solCont = NDSolveProcessSolutions[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]

• What a well founded answer! Thanks. Where did you find these interesting variants of NDSolve? Some weeks ago I tried to use NDSolve in a given interpolation grid to speed up simulation time. Do you think it's possible with NDSolveProcessEquations to realize something like "InterpolationPointsSubdivision"? – Ulrich Neumann Dec 29 '20 at 21:22
• @UlrichNeumann Some years ago I decided I wanted to be familiar with NDSolve, so I downloaded Advanced Numerical Differential Equation Solving in the Wolfram Language and read through it like it was a juicy mystery novel. -- With something like stops = Union@Flatten@Cases[ode, if_InterpolatingFunction :> if@"Grid, Infinity] and the Do loop above, one should be able to implement an "InterpolationPointsSubdivision" method. – Michael E2 Dec 29 '20 at 22:28
• Thanks! This will help a lot. It turns out that I have found another straightforward way to resolve the problem: to divide both the jacobian and rate functions by the same large value, say, 10^24: rate[T_] = 10^Interpolation[datatable[[All, {1, 2}]], InterpolationOrder -> 1][ T]/10^24; jacobian[T_] = Interpolation[datatable[[All, {1, 3}]], InterpolationOrder -> 1][T]/ 10^24; However, the approach you have proposed should work for more complicated systems, which I have already faced. – John Taylor Dec 29 '20 at 23:29
• @MichealE2 Thank you very much, I' ll download this novel too and try your suggestion! – Ulrich Neumann Dec 30 '20 at 10:32
• @JohnTaylor Your straightforward way "scales" your ode in the right way. Often this scaling might be founded due to physical/mathematical reasons(In this case for example change of time scaling ) – Ulrich Neumann Dec 30 '20 at 10:37

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.