NDSolve with varying PrecisionGoal and WorkingPrecision

Sometimes we need higher numerical precision to deal with large number cancellation in an equation. But if this cancellation happens only in a small (and known) parameter space, would it be possible to only use high precision calculation in that small part of the parameter space? For example, something like

NDSolve[{y'[x] == y[x] Cos[x + y[x]], y == 1}, y, {x, 0, 30},
WorkingPrecision -> If[x < 5, $MachinePrecision, 50], PrecisionGoal -> If[x < 5,$MachinePrecision, 50]]

This code doesn't work, because Mathematica doesn't relate the x in the two If expressions to the x in the ODE. Is there any way to make the above idea work?

• Note that $MachinePrecision is different from MachinePrecision. WorkingPrecision ->$MachinePrecision uses arbitrary precision with $MachinePrecision digits of precision, which may be what you wanted. One can use WorkingPrecision -> MachinePrecision, which is the default, for machine precision calculation, when desired. – Michael E2 Mar 2 '14 at 13:52 1 Answer No direct way I'm aware of to do this: as far as I know, when the initial NDSolve is evaluated, those options a plugged into an internal slot and can't be tinkered. However, perhaps something like this would serve your purpose: low = NDSolve[{y'[x] == y[x] Cos[x + y[x]], y == 1}, y, {x, 0, 5}, WorkingPrecision ->$MachinePrecision,
PrecisionGoal -> \$MachinePrecision];

high = NDSolve[{y'[x] == y[x] Cos[x + y[x]], y == 1},
y, {x, 5, 30}, WorkingPrecision -> 50, PrecisionGoal -> 50];

myinterpolating[x_] =
Piecewise[{{(y /. First@low)[x], x < 5}, {(y /. First@high)[x],
x >= 5}}];

Column[{x = myinterpolating, Precision[x], y = myinterpolating,
Precision[y]}]

(*

0.528357862862939
15.3332
0.06434904783198493788519509130059592424856983961826
49.0072

*)
• Thanks for the answer! One thing I was not sure is, for high, because there is y == 1 condition, does NDSolve really start from x=5, or it has to start from x=0 and thus has the same performance as starting to integrate using high performance from the beginning? – Yi Wang Mar 2 '14 at 9:23
• @YiWang: Good question - I can't answer definitively (code-speak for "beats me"). I suppose some timings are in order...and I'm sure experts in the function will chime in. Edit: just did some quick timings, it takes longer with a higher lower bound - puzzling, so I'm not sure if my idea will be of any use to you. – ciao Mar 2 '14 at 9:27
• I think this doesn't really do what you expect and the OP intended: if you give a start value for x=0 but ask NDSolve to solve for x=5 to x=30 in the second call it will actually solve from x=0 with high precision but only store the results from 5 to 30 in the result. If you want to only solve with low precision you should instead give the lower precision endpoint result as starting value for the second call, something like: y == SetPrecision[(y /. First@low), 50]... – Albert Retey Mar 2 '14 at 10:53
• @AlbertRetey: I thought of that but didn't try it - working on other stuff. Still puzzling that it takes more time with a higher lower bound than just doing the whole range. Perhaps combining that code with mine will be useful to OP, I'll tinker with the idea when time permits. – ciao Mar 2 '14 at 11:08
• I'm not sure but I think that the difference might be that with a higher lower bound NDSolve will try to make an exact match in the x-steps for that value which might result in a slightly different discretization. You can see similar effects with e.g. NDSolve[{y'[x] == y[x] Cos[x + y[x]], y == 1}, y, {x, 0, 5, 30}, WorkingPrecision -> 50, PrecisionGoal -> 50]. See e.g. this answer about the undocumented feature I used here for the x range specification... – Albert Retey Mar 2 '14 at 13:56