Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
    
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 at 13:52

1 Answer 1

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

*)
share|improve this answer
    
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? –  Yi Wang Mar 2 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. –  rasher Mar 2 at 9:27
1  
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]... –  Albert Retey Mar 2 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. –  rasher Mar 2 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[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... –  Albert Retey Mar 2 at 13:56

Your Answer

 
discard

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.