2
$\begingroup$

I am trying to solve a differential equation that contains some interpolating functions as input, and for reasons unknown to me I am encountering an error which makes Mathematica use extrapolation. I managed to create a toy problem which seems to reproduce my error below, and I hope someone can explain me why I am getting the error I am getting.

The following numerically solves a simple differential equation:

   Asol = NDSolve[{A'[r] + A[r] == 0, A[1] == 1}, A, {r, 1, 0}, WorkingPrecision -> 30]

and gives me an interpolating function as output:

   {{A -> InterpolatingFunction[{{0, 1.00000000000000000000000000000}}, <>]}}.

Now I would like to numerically solve a differential equation that contains this interpolating function:

   NDSolve[{(B'[r] + Sqrt[A[r]] /. Asol) == 0, B[1] == 2}, B, {r, 0, 1}].

Although Mathemathica manages to solve it, it does so while giving the following error:

  InterpolatingFunction::dmval: Input value {1.} lies outside the range 
  of data in the interpolating function. Extrapolation will be used. >>

  {{B -> InterpolatingFunction[{{0., 1.}}, <>]}}

When I copy the input value in the error (i.e. the {1.}) and paste it, it reveals that the value is actually 1.0000000149011612`, which is indeed not in the range of the interpolating function.

My question is: why does Mathematica try to evaluate it at this point in the first place, as it is outside the range I am requesting in NDSolve? And how can I prevent it?

Some comments/things I noticed:

  • The error does not seem to occur if I try to solve the equation B'[r]+A[r]==0. It only occurs when I the equation contains the square root, or any non-integer power of A[r].
  • The error also does not seem to occur if I remove the WorkingPrecision -> 30 part in the command Asol. However, it does occur even if I also add WorkingPrecision -> 30 to the differential equation for B.
  • Of course I can solve these equations exactly, but in my actual program I am using analogues of A and B which do not have exact solutions. In this program I would like a quite precise answer, therefore I do need the WorkingPrecision command.
$\endgroup$
  • $\begingroup$ I would guess that NDSolve is going beyond the upper bound of integration due to round-off issues of some sort. Of course, only individuals familiar with the guts of NDSolve would be able to answer with certainty. In any case, this should have no impact on the accuracy of the final answer. Turn off this message, if it bothers you. $\endgroup$ – bbgodfrey Jul 18 '17 at 2:24
  • $\begingroup$ @bbgodfrey Thanks for your comment. My only worry is that in my actual programme it sometimes happens that the differential equation is too hard to solve, and then this error pop ups as well, in which case it is kind of serious. That's why I wouldn't turn it off. Anyway, I'll just try to solve for A in a bigger range than for B, that might help. $\endgroup$ – ScroogeMcDuck Jul 18 '17 at 8:38
  • $\begingroup$ My guess is that it's taking a step (or more) forward to approximate the derivatives for the sake of error estimation. If so, there's probably not a lot you can do conveniently with your setup to avoid it. (It could step backwards, but I don't know of an option or other trick that will let you do that.) $\endgroup$ – Michael E2 Mar 15 at 21:29
  • $\begingroup$ Some evidence: if you reflect r -> -u, then the warning goes away: NDSolve[{(-BB'[u] + Sqrt[A[-u]] /. Asol) == 0, BB[-1] == 2}, BB, {u, -1, 0}] -- BTW, the warning probably does not matter, since it the steps outside the interval are so small. Another BTW, you will get a more accurate interpolating function for A, if you use the option InterpolationOrder -> All in NDSolve. $\endgroup$ – Michael E2 Mar 15 at 21:44
1
$\begingroup$

Instead of doing "nested" NDSolve applications, include both equations in a single NDSolveValue:

sol = NDSolveValue[
    {
    A'[r] + A[r] == 0, A[1] == 1,
    B'[r] + Sqrt[A[r]] == 0, B[1] == 2
    },
    {A, B},
    {r, 0, 1},
    WorkingPrecision->30
];

Visualization:

Plot[Evaluate @ Through @ sol[t], {t, 0, 1}, PlotLabels -> {A, B}]

enter image description here

$\endgroup$
  • $\begingroup$ I'd like to add that this is the preferred way to do this, because nesting NDSolve will propagate errors forwards. It's better to solve the whole thing in one go to make sure Mathematica can track the convergence properly. $\endgroup$ – Sjoerd Smit Mar 15 at 21:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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