3
$\begingroup$

I was considering the following piece of code

u[y_, z_] = 36*y*(1 - y)*z*(1 - z);
sys = {x'[t] == u[y[t], z[t]], y'[t] == 0, z'[t] == -2.2*^-3, 
   x[0] == 0, y[0] == y0, z[0] == z0, 
   WhenEvent[x[t] > 500, "StopIntegration"]};
sol = Table[
   NDSolve[sys, {x, y, z}, {t, 0, 500}], {y0, 0.1, 0.9, 0.1}, {z0, 
    0.1, 0.9, 0.1}];
ParametricPlot3D[{x[t], y[t], z[t]} /. sol, {t, 0, 500}, 
 PlotRange -> {{0, 1000}, {0, 1}, {0, 1}}, ImageSize -> 500, 
 BoxRatios -> {5, 1, 1}]

However, it does not seem that WhenEvent captures correctly the event x>500, indeed Mathematica does not stop integration. I've also tried with "==" but I have the same problem.

What is wrong?

Thanks!

Improve this question
$\endgroup$
4
  • 2
    $\begingroup$ I think it's actually working. You're just plotting the solutions for t from 0 to 500 even for the solutions that stop at some earlier time. Try looking at just the solution for {y0 -> 0.5 , z0 -> 0.9} and do x["Domain"]/.sol to get the domain of the function. Also, you can replace "StopIntegration" with Print[{y0, z0, t}]; "StopIntegration" to make Mathematica spit out the times at which it stops and for which parameters. $\endgroup$
    – march
    Commented Oct 22, 2015 at 16:02
  • 1
    $\begingroup$ The core issue here is Plot (et al) suppress the warnings about extrapolation of the InterploatingFunction. Is there a way to disable that (ie. show the warnings )? $\endgroup$
    – george2079
    Commented Oct 22, 2015 at 16:56
  • $\begingroup$ @March: Yes, I've realized just now that it is a problem of extrapolation outside the integration interval. However, Print[{y0,z0,t}] does not print anything, furthermore x["Domain"]/.sol gives {0,500} and it seems to me that the integration does not stop at all. SORRY! I was unaware that when you put in more events these must be listed! Actually it is working, thanks! $\endgroup$
    – Fabio
    Commented Oct 22, 2015 at 18:38
  • $\begingroup$ Any clue how to Reap/Sow in WhenEvent? $\endgroup$
    – Fabio
    Commented Oct 22, 2015 at 18:45

1 Answer 1

Reset to default
3
$\begingroup$ 
u[y_, z_] := 36*y*(1 - y)*z*(1 - z);
sys = {
   x'[t] == u[y[t], z[t]], y'[t] == 0,
   z'[t] == -2.2*^-3,
   x[0] == 0, y[0] == y0, z[0] == z0, 
   WhenEvent[x[t] == 500, "StopIntegration"]};

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
sol = Table[NDSolveValue[sys, {x, y, z}, {t, 0, 500}], {y0, 0.1, 0.9, 0.2}, 
                                                       {z0, 0.1, 0.9, 0.2}];
doms = InterpolatingFunctionDomain[#[[1]]][[1, 2]] & /@  Flatten[sol, 1];
tt = Transpose[{Flatten[sol, 1], doms}];
plots = ParametricPlot3D[Through[#[[1]][t]], {t, 0, #[[2]]}, 
                         BoxRatios -> {5, 1, 1}] & /@ tt;
prs = PlotRange /. AbsoluteOptions[#, PlotRange] & /@ plots;
prs1 = {Min[#[[1]]], Max[#[[2]]]} & /@ Transpose /@ Transpose@prs;
Show[plots, PlotRange -> prs1]

Mathematica graphics

edit

You may plot the ending time wrt y0 and z0:

coords = Table[{y0, z0}, {y0, 0.1, 0.9, 0.1}, {z0, 0.1, 0.9, 0.1}];
prs = (Last@First@(PlotRange /. AbsoluteOptions[#, PlotRange])) & /@  plots;
ListPlot3D[Join @@@ Transpose@{Flatten[coords, 1], List /@ prs}]

Mathematica graphics

Improve this answer
$\endgroup$
1
  • $\begingroup$ Interesting solution though complicate. I have to study in detail, thanks. $\endgroup$
    – Fabio
    Commented Oct 22, 2015 at 18:35

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.