3
$\begingroup$

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[0] == 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?

Improve this question
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – Michael E2
    Commented Mar 2, 2014 at 13:52

1 Answer 1

Reset to default
3
$\begingroup$

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[0] == 1}, y, {x, 0, 5}, 
   WorkingPrecision -> $MachinePrecision, 
       PrecisionGoal -> $MachinePrecision];

high = NDSolve[{y'[x] == y[x] Cos[x + y[x]], y[0] == 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[2], Precision[x], y = myinterpolating[10],
   Precision[y]}]

(*

0.528357862862939 
15.3332
0.06434904783198493788519509130059592424856983961826
49.0072

*)
Improve this answer
$\endgroup$
5
  • $\begingroup$ Thanks for the answer! One thing I was not sure is, for high, because there is y[0] == 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? $\endgroup$
    – Yi Wang
    Commented Mar 2, 2014 at 9:23
  • $\begingroup$ @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. $\endgroup$
    – ciao
    Commented Mar 2, 2014 at 9:27
  • 1
    $\begingroup$ 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[0] == SetPrecision[(y /. First@low)[5], 50]... $\endgroup$
    – Albert Retey
    Commented Mar 2, 2014 at 10:53
  • $\begingroup$ @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. $\endgroup$
    – ciao
    Commented Mar 2, 2014 at 11:08
  • $\begingroup$ 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[0] == 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... $\endgroup$
    – Albert Retey
    Commented Mar 2, 2014 at 13:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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