1
$\begingroup$

In Mathematica 9.0.1, here is a simple code to numerically solve an initial-value problem:

sol = With[{k = 0.4, x0 = .01, xMax = 20}, 
  ParametricNDSolve[
    {y''[x] + 3/x y'[x] - y[x] + 3/2 y[x]^2 - k/2 y[x]^3 == 0 (*diff. eq.*), 
     y[x0] == y0 + 1/8 x0^2 (y0 - 3/2 y0^2 + k/2 y0^3),       (*init. cond. 1*)
     y'[x0] == 1/4 x0 (y0 - 3/2 y0^2 + k/2 y0^3)              (*init. cond. 2*), 
    WhenEvent[y[x] == 0, "StopIntegration"]}, y, {x, x0, xMax}, y0]
  ]

Then I can make a plot of a solution (with y0=4.4), and I get:

Plot[y[4.4][x] /. sol, {x, .01, 20}, PlotRange -> {-1, 5}]

enter image description here

But if I modify the WhenEvent in the code above so that it reads:

WhenEvent[y'[x] == 0, "StopIntegration"]}, y, {x, x0, xMax}, y0],

the solution to the differential equation becomes different:

enter image description here

I expected NDSolve to stop integrating when the curve begins to turn around near x=5. Instead, the whole solution changes -- and it doesn't stop integrating. Am I misinterpreting WhenEvent?

$\endgroup$
  • 2
    $\begingroup$ I'm not sure what you are doing, but for me y[4.4]/.sol returns an InterpolatingFunction which is valid from 0.01 to 5.14 and that seems to be the value where y' becomes zero as desired. You can of course still plot it but for values larger than x=5.14 but then only an extrapolation of that InterpolatingFunction is plotted, nothing that has anything to do with the real solution. I'm not sure why they changed the older behaviour where you'd have seen a bunch of warning messages but in any case you should check the end point before plotting when using "StopIntegration"... $\endgroup$ – Albert Retey May 20 '15 at 19:35
  • $\begingroup$ Headslap! Yes.. you're right. $\endgroup$ – QuantumDot May 20 '15 at 19:37
  • $\begingroup$ great, I wasn't sure but as I struggled about this myself it was an easy guess :-) $\endgroup$ – Albert Retey May 20 '15 at 19:39
  • $\begingroup$ I have added my comment as an answer, so if you want you can accept it and make that question disappear from the unanswered list... $\endgroup$ – Albert Retey May 21 '15 at 7:44
0
$\begingroup$

If you look at the output of:

y[4.4]/.sol

you will find that that is an interpolating function which is valid from x=0.01 to roughly 5.14. If one plots that result from x=0.01 to 5.14 that shows that 5.14 is the point where the derivative becomes zero. So it looks like NDSolve does the right thing. What happens when you Plot from x=0.01 to xMax=20 is that Mathematica uses an extrapolation of the result for values larger than 5.14 in the plot, and of course that extrapolation doesn't have to anything with the real solution for values larger than 5.14. Older versions of Mathematica would have given some warning messages about having used extrapolations of an interpolating function, but newer don't do that anymore. I'm not sure whether that was a wise design decision as these messages would have helped you to immediately see the problem and in the probably relatively unlikely case you realy wanted to use extrapolations it would be easy enough to switch those messages off.

Anyway, if ever you use a WhenEvent with "StopIntegration" you better extract the integration end point and use that when making use of the result, e.g. for a plot:

With[{xMax = ((y[4.4] /. sol)@"Domain")[[1, 2]]},
  Plot[y[4.4][x] /. sol, {x, .01, xMax}, PlotRange -> {-1, 5}]
]

I would also strongly recommend to learn about the new version of *NDSolveValue call functions which directly return the interpolating functions, which makes such operations somewhat more convenient:

sol = With[{k = 0.4, x0 = .01, xMax = 20}, 
  ParametricNDSolveValue[{
     y''[x] + 3/x y'[x] - y[x] + 3/2 y[x]^2 - k/2 y[x]^3 == 0 (*diff.eq.*), 
     y[x0] == y0 + 1/8 x0^2 (y0 - 3/2 y0^2 + k/2 y0^3),    (*init.cond.1*)
     y'[x0] ==  1/4 x0 (y0 - 3/2 y0^2 + k/2 y0^3)          (*init.cond.2*), 
     WhenEvent[y'[x] == 0, "StopIntegration"]
   }, y, {x, x0, xMax}, y0]
]

 Plot[sol[4.4][x], {x, .01, (sol[4.4]@"Domain")[[1, 2]]}, 
   PlotRange -> {-1, 5}]
$\endgroup$

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.