Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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] /. %


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?

share|improve this question
up vote 6 down vote accepted

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]
share|improve this answer
Thanks a lot! This indeed fixed the problem. :D – stupidity May 26 '14 at 2:21
@stupidity You're welcome! – Michael E2 May 26 '14 at 2:22

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.