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.

I have a system of ODEs that I solve. During the integration process, there's an event that I want to catch, but I want to (a) continue the integration after the event and (b) catch only the first one.

For example, let's say I want to calculate $\pi/2$ by this method:

flag = False;
NDSolve[{y''[t] == -y[t], y[0] == 0, y'[0] == 1}, y, {t, 0, 20}, 
 Method -> {"EventLocator", "Event" -> ((! flag) && y'[t] == 0), 
   "EventAction" :> (flag = True; piHalf = t;)}]

For some reason this does not work - the events keep being triggered even though the value of flag does change.

screenshot

In the end of the process I get $\pi/2=17.28$ (this might be disastrous if I'll use this value to build airplanes!) How can I solve this?

share|improve this question
add comment

2 Answers

up vote 6 down vote accepted

I think there are two misunderstandings here, both are somewhat subtle. First, giving the "Event" option as you did doesn't work because it is immediately evaluated, so that NDSolve does only see this:

flag = False;
And[! flag, y'[t] == 0]

(*
==> Derivative[1][y][t] == 0
*)

meaning it never even sees that the event should depend on flag at all. This can be circumvented by using RuleDelayed (:>) instead of Rule (->). Unfortunately the result will still not be what you expect. The reason is that y'[t]==0 will only fire if at any step y'[t]==0 evaluates to True, which will only happen by chance and/or if the step size is small enough. It is much better/more reliable to let NDSolve check for changes in sign and the combination of these changes will make your example work as I think you expected:

flag = False; piHalf =.;
NDSolve[{y''[t] == -y[t], y[0] == 0, y'[0] == 1}, y, {t, 0, 20}, 
  Method -> {
    "EventLocator",
    "Event" :> And[! flag, y'[t] >= 0], 
    "EventAction" :> (flag = True; piHalf = t;)
    }
  ];
piHalf

I just have seen Nassers answer and wanted to add some remarks about it:

1) Nassers WhenEvent basically is equivalent to these "EventLocator" settings:

flag = False; piHalf =.;
NDSolve[{y''[t] == -y[t], y[0] == 0, y'[0] == 1}, y, {t, 0, 20}, 
  Method -> {
    "EventLocator",
    "Event" -> y'[t],
    "EventAction" :> (If[! flag, piHalf = t]; flag = True;)
    }
  ];
piHalf

which is somewhat different: the "EventAction" is in this case evaluated at every sign change of y'[t], but only the first time piHalf is set. The other solution will only evaluate "EventAction" once. For the example at hand the net result (value of piHalf) is of course the same.

2) WhenEvent has the attribute HoldAll, which makes it behave more like the RuleDelayed case of the "Event" setting which solves the evaluation problem. Using a predicate in WhenEvent which is just refering to a global variable seems to cause kernel crashes, but there is a very elegant way to achieve that the event will only be triggered once without even bothering with a global variable:

NDSolve[{y''[t] == -y[t], y[0] == 0, y'[0] == 1,
   WhenEvent[y'[t] == 0, piHalf = t; "RemoveEvent"]}, y, {t, 0, 20}
  ];
piHalf

3) if you decide to go with the "EventLocator" method you might also want to have a look at the option "Direction" which lets you discriminate between sign changes from below vs. from above.

share|improve this answer
    
@Nasser: this was not meant as critisism of your solution, the original question wasn't very percise about which of the two would be prefered, if any. I just wanted to explain in how the two solutions differ and also to justify that I wrote another answer. In fact your solution made me look closer into how to use WhenEvent and caused me to learn about "RemoveEvent" myself, so thanks for that :-). –  Albert Retey Dec 18 '12 at 17:02
add comment

I do not know why EventAction is not working as you show. (could be a bug?) but I just want to give a solution using WhenEvent (new in V9) which does the same thing you are trying to do which is to assign a value to the variable piHalf only once (first time the event happens), but continue the NDSolve action.

flag = False;
sol = NDSolve[{y''[t] == -y[t], y[0] == 0, y'[0] == 1,
   WhenEvent[y'[t] == 0,
    If[Not[flag], flag = True; piHalf = t; Print[flag, " ", " t=", t]]
    ]
   },
  y, {t, 0, 20}
  ]

(* True   t=1.5707963620260097
   Out[4]= {{y -> InterpolatingFunction[]}} *)

piHalf
(*1.5708*)

Plot[Evaluate[y[t] /. First[sol]], {t, 0, 20},
     Frame -> True, 
     FrameLabel -> {{y[t], None}, {t, "capturing event first time"}},
     GridLines -> Automatic,
     Epilog -> {Red, PointSize[Large], Point[{piHalf, 0}]}
    ]

Here is a screen shot of the whole process

Mathematica graphics

share|improve this answer
    
Thanks @Nasser. Unforutnately, though, I don't have V9, so i can't use that solution. I found an ad-hoc workaround that I'll post soon, unless someone else gives a good solution. –  yohbs Dec 18 '12 at 16:18
add comment

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.