differential equations - NDSolve::ndsz problem

I'm using NDSolve to solve a set of coupled differential equation which depend on a variable x. I noticed that when I set the range of x from a small value to a large value, I obtain solutions. But when I do the opposite I get a problem!

For example when I consider this,

t1 = 3000;
t2 = 4*^16;
f = 1/(16 Pi^2);
NDSolve[{y'[x] == f/x * 16 y[x]^3, y[t1] == 0.37}, y, {x, t1, t2},
Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4}];

No errors and I can find:

y[t2] /. %


{0.920385}

Then I use this value of y[t2] as the initial condition for the following case, where I swap the range of x:

t1 = 3000;
t2 = 4*^16;
f = 1/(16 Pi^2);
NDSolve[{y'[x] == f/x * 16 y[x]^3, y[t2] == 0.92}, y, {x, t2, t1},
Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4}]

I get the error:

NDSolve::ndsz: At x == 4.`*^16, step size 
is effectively zero; singularity  or stiff system suspected.

I tried to switch off the StiffnessTest to see if the problem goes away, which would mean that it's a stiffness problem, but the error message still shows. Also tried to use:

Method -> {StiffnessSwitching, Method -> {ExplicitRungeKutta, Automatic}}

And it didn't work either.

Considering that this is a typical differential equation encountered in physics (renormalization group equations) I'm pretty sure that there should be no singularity at x = t2, and I should be able to solve from the high-scale to the low-scale.

Any insights on why this is happening and how to deal with it?

Answer

Probably it's round-off error. Increase WorkingPrecision, and limit AccuracyGoal and PrecisionGoal. Just what to set them to depends on the actual equations.

This works for the OP's example:

t1 = 3000;
t2 = 4*^16;
f = 1/(16 Pi^2);
sol1 = NDSolve[{y'[x] == f/x*16 y[x]^3, y[t1] == 0.37`50}, 
   y, {x, t1, t2}, 
   Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4}, 
   WorkingPrecision -> 50, MaxSteps -> 20000, AccuracyGoal -> 15, 

   PrecisionGoal -> 15];

sol2 = NDSolve[{y'[x] == f/x*16 y[x]^3, 
  y[t2] == (y[t2] /. First@sol1)}, y, {x, t1, t2}, 
 Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4}, 
 WorkingPrecision -> 50, AccuracyGoal -> 15, PrecisionGoal -> 15, 
 MaxSteps -> 20000]