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.
Improve this question
$\endgroup$
4
  • $\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, 2017 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, 2017 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, 2019 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, 2019 at 21:44

1 Answer 1

Reset to default
3
$\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

Improve this answer
$\endgroup$
1
  • $\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, 2019 at 21:21

Your Answer

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