|
EventLocator IntroductionIt is often useful to be able to detect and precisely locate a change in a differential system. For example, with the detection of a singularity or state change, the appropriate action can be taken, such as restarting the integration. An "event" for a differential system  is a point along the solution at which a real valued event function  is zero. It is also possible to consider boolean valued event functions, in which case the event occurs when the function changes from True to False or vice-versa. The EventLocator method which is built into NDSolve works effectively as a controller method; it handles checking for events and taking the appropriate action, but the integration of the differential system is otherwise left completely to an underlying method. In this section, examples are given to demonstrate the basic use of the EventLocator method and options. Subsequent sections show more involved applications of event location, such as period detection, Poincaré sections, and discontinuity handling. These initialization commands load some useful packages that have some differential equations to solve and define some utility functions. In[6]:=  |
A simple example is locating an event, such as when a pendulum started at a non-equilibrium position will swing through its lowest point and stopping the integration at that point. This integrates the pendulum equation up to the first point at which the solution x[t] crosses the axis. In[5]:=  |
Out[5]=
|
From the solution you can see that pendulum reaches its lowest point x[t] = 0 at about t = 1.675. Using the InterpolatingFunctionAnatomy package, it is possible to extract the value from the InterpolatingFunction object. This extracts the point at which the event occurs and makes a plot of the solution (black) and its derivative (blue) up to that point. In[6]:=  |
Out[7]=
|
When you use the event locator method, the events to be located and the action to take upon finding an event are specified through method options of the EventLocator method. The default action on detecting an event is to stop the integration as demonstrated above. The event action can be any expression. It is evaluated with numerical values substituted for the problem variables whenever an event is detected. This prints the time and values each time the event x'[t] = x[t] is detected for a damped pendulum. In[8]:=  |
Out[8]=
|
Note that in the above example, the "EventAction" option was given using RuleDelayed (⧴) to prevent it from evaluating except when the event is located. You can see from the printed output that when the action does not stop the integration, multiple instances of an event can be detected. Events are detected when the sign of the event expression changes. You can restrict the event to be only for a sign change in a particular direction using the "Direction" option. This collects the points at which the velocity changes from negative to positive for a damped driven pendulum. Reap and Sow are programming constructs which are useful for collecting data when you don't, at first, know how much data there will be. Reap[expr] gives the value of expr together with all expressions to which Sow has been applied during its evaluation. Here Reap encloses the use of NDSolve and Sow is a part of the event action, which allows us to collect data for each instance of an event. In[9]:=  |
Out[9]=
|
You may notice from the output of the previous example that the events are detected when the derivative is only approximately zero. When the method detects the presence of an event in a step of the underlying integrator (by a sign change of the event expression), then it uses a numerical method to approximately find the position of the root. Since the location process is numerical, you should expect only approximate results. Location method options AccuracyGoal, PrecisionGoal, and MaxIterations can be given to those location methods which use FindRoot to control tolerances for finding the root. It is also desirable to be able to detect events which are the result of some modification which is external to the differential system itself, such as a user pressing a button. It is often hard, if not impossible to represent such events as continuous functions with roots at the correct point. Often, it is much more natural to reflect such a state change by a function which is either True or False. For boolean valued event functions, an "event" occurs when the function switches from True to False or vice-versa. The "Direction" option can be used to restrict the event only from changes from True to False ("Direction"->-1) or only from changes from False to True ("Direction"->1). This opens up a small window with a button, which when pressed changes the value of the variable stop to True from its initialized value of False. In[8]:=  |
This integrates the pendulum equation up until the button is pushed (or the system runs out of memory) In[6]:=  |
Out[6]=
|
Take note that in this example, the "Event" option was specified with RuleDelayed (:>) to prevent the immediate value of stop from being evaluated and set up as the function. You can specify more than one event. If the event function evaluates numerically to a list, then each component of the list is considered to be a separate event. You can specify different actions, directions, etc. for each of these events by specifying the values of these options as lists of the appropriate length. This integrates the pendulum equation up until the point at which the above "Stop" button is pressed. The number of complete swings of the pendulum is kept track of during the integration. In[7]:=  |
Out[8]=
|
As you can see from the previous example, it is possible to mix real and boolean valued event functions. The expected number of components and type of each component is based on the values at the initial condition and needs to be consistent throughout the integration. The Method option of EventLocator allows the specification of the numerical method to use in the integration. Event Location MethodsThe EventLocator method works by taking a step of the underlying method and checking to see if the sign (or parity) of any of the event functions is different at the step endpoints. Event functions are expected to be real or boolean valued, so if there is a change, there must be an event in the step interval. For each event function which has an event occurrence in a step, a refinement procedure is carried out to locate the position of the event within the interval. There are several different methods which can be used refine the position. These include simply taking the solution at the beginning or the end of the integration interval, a linear interpolation of the event value and using bracketed root finding methods. The appropriate method to use depends on a trade off between execution speed and location accuracy. If the event action is to stop the integration then the particular value at which the integration is stopped depends on the value obtained from the EventLocationMethod option of EventLocator. Location of a single event is usually fast enough so that the method used will not significantly influence the overall computation time. However, when an event is detected multiple times, the location refinement method can have a substantial effect. StepBegin and StepEnd MethodsThe crudest methods are appropriate for when the exact position of the event location doesn't really matter or does not reflect anything with precision in the underlying calculation. The stop button example from the previous section is such a case: time steps are computed so quickly that there is no way that one can time the press of a button to be within a particular time step, much less at a particular point within a time step. Thus, based on the inherent accuracy of the event, there is no point in refining at all. You can specify this by using the "StepBegin" or "StepEnd" location methods. In any example where the definition of the event is heuristic or somewhat imprecise, this can be an appropriate choice. LinearInterpolation MethodWhen event results are needed for the purpose of points to plot in a graph, you only need to locate the event to the resolution of the graph. While just using the step end is usually too crude for this, a single linear interpolation based on the event function values suffices. You can specify refinement based on a single linear interpolation with the setting "LinearInterpolation". This computes the solution for a single period of the pendulum equation and plots the solution for that period. In[9]:=  |
Out[10]=
|
At the resolution of the plot over the entire period, you cannot see that the endpoint may not be exactly where the derivative hits the axis. However, if you zoom in enough, you can see the error. This shows a plot just near the endpoint. In[11]:=  |
Out[11]=
|
The linear interpolation method is sufficient for most viewing purposes, such as the Poincaré section examples shown in the following section. Note that for Boolean valued event functions, linear interpolation is effectively only one bisection step, so the linear interpolation method may be inadequate for graphics. Brent's MethodThe default location method is the event location method "Brent", finding the location of the event using FindRoot with Brent 's method. Brent's method starts with a bracketed root and combines steps based on interpolation and bisection, guaranteeing a convergence rate at least as good as bisection. You can control the accuracy and precision to which FindRoot tries to get the root of the event function using method options for the "Brent" event location method. The default is to find the root to the same accuracy and precision as NDSolve is using for local error control. For methods which support continuous or dense output, the argument for the event function can be found quite efficiently simply by using the continuous output formula. However, for methods which do not support continuous output, the solution needs to be computed by taking a step of the underlying method, which can be relatively expensive. An alternate way of getting a solution approximation which is not accurate to the method order, but is consistent with using FindRoot on the InterpolatingFunction object returned from NDSolve is to use cubic Hermite interpolation, obtaining approximate solution values in the middle of the step by interpolation based on the solution values and solution derivative values at the step ends. ComparisonThis example integrates the pendulum equation for a number of different event location methods and compares the time when the event is found. This defines the event location methods to use. This integrates the system and prints out the method used and the value of the independent variable when the integration is terminated.
 ExamplesFalling BodyThis system models a body falling under the force of gravity encountering air resistance (see [M04]). The event action stores the time when the falling body hits the ground and stops the integration. In[11]:=  |
Out[11]=
|
This plots the solution and highlights the initial and final points (green and red) by encircling them. In[15]:=  |
Out[17]=
|
Period of the Van der Pol OscillatorThe Van der Pol oscillator is an example of an extremely stiff system of ODEs. The event locator method can call any method for actually doing the integration of the ODE system. The default method, Automatic, automatically switches to a method appropriate for stiff systems when necessary, so that stiffness does not present a problem. This integrates the Van der Pol system for a particular value of the parameter  = 1000 up to the point where the variable  reaches its initial value and direction. In[20]:=  |
Out[20]=
|
Note that the event at the initial condition is not considered. By selecting the endpoint of the NDSolve solution, it is possible to write a function which returns the period as a function of . This defines a function which returns the period as a function of . In[21]:=  |
This uses the function to compute the period at = 1000 In[22]:=  |
Out[22]=
|
Of course, it is easy to generalize this method to any system with periodic solutions. Poincaré SectionsThe Hénon-Heiles systemThis is a useful technique for visualizing the solutions of high dimensional differential systems. Define the Hénon-Heiles system which models stellar motion in a galaxy. This gets the Hénon-Heiles system from the NDSolveProblems package. In[23]:=  |
Out[24]=
|
The Poincaré section of interest in this case is the collection of points in the  plane when the orbit passes through  . Since the actual result of the numerical integration is not required, it is possible to avoid storing all the data in InterpolatingFunctions by specifying the output variables list (in the second argument to NDSolve) as empty, or {}. This means that NDSolve will produce no InterpolatingFunctions as output, avoiding storing a lot of unnecessary data. NDSolve does give a message NDSolve::noout warning there will be no output functions, but it can safely be turned off in this case since the data of interest is collected from the event actions. The linear interpolation event location method is used because the purpose of the computation here is to view the results in a graph with relatively low resolution. If you were doing an example where you needed to zoom in on the graph to great detail or find a feature, such as a fixed point of the Poincaré map, it would be more appropriate to use the default location method. This turns off the message warning about no output. In[25]:=  |
This integrates the Hénon-Heiles system using a fourth order explicit Runge Kutta method with fixed step size of 0.25. The event action is to Sow the values of  In[26]:=  |
This plots the Poincaré section. The collected data is found in the last part of the result of Reap and the list of points is the first part of that. In[31]:=  |
Since the Hénon-Heiles system is Hamiltonian, a symplectic method gives much better qualitative results for this example. This integrates the Hénon-Heiles system using a fourth order symplectic partitioned Runge Kutta method with fixed step size of 0.25. The event action is to Sow the values of  In[33]:=  |
This plots the Poincaré section. The collected data is found in the last part of the result of Reap and the list of points is the first part of that. In[34]:=  |
The ABC flowThis loads an example problem of the Arnold-Beltrami-Childress (ABC) flow which is used to model chaos in laminar flows of the three-dimensional Euler equations. In[25]:=  |
This defines a splitting Y1, Y2 of the system by setting some of the right-hand side components to zero. In[29]:=  |
Out[29]=
|
In[30]:=  |
Out[30]=
|
This defines the implicit midpoint method. In[31]:=  |
This constructs a second-order splitting method that retains volume and reversing symmetries. In[32]:=  |
This defines a function which gives the Poincaré section for a particular initial condition. In[34]:=  |
This finds the Poincaré sections for several different initial conditions and flattens them together into a single list of points. In[35]:=  |
Bouncing BallThis example is a generalization of an example in [SGT03]. It models a ball bounding down a ramp with a given profile. The example is good for demonstrating how you can use multiple invocation of NDSolve with event location to model some behavior. This defines a function which computes the solution from one bounce to the next. The solution is computed up until the next time the path intersects the ramp. In[36]:=  |
This defines the function for the bounce when the ball hits the ramp. The formula is based on reflection about the normal to the ramp assuming only the fraction k of energy is left after a bounce. In[37]:=  |
This defines the function which runs the bouncing ball simulation for a given reflection ratio, ramp, and starting position. In[38]:=  |
This is the example which is done in [SGT03]. In[39]:=  |
The ramp is now defined to be a quarter circle In[41]:=  |
This adds a slight waviness to the ramp. In[43]:=  |
Event DirectionThis example illustrates the solution of the restricted two body problem, a standard nonstiff test system of four equations. The example traces the path of a spaceship travelling around the moon and returning to the earth (see p. 246 of [SG75]). The ability to specify multiple events and the direction of the zero crossing is important. The initial conditions have been chosen to make the orbit periodic. The value of corresponds to a spaceship travelling around the moon and the earth. The event function is the derivative of the distance from the initial conditions. A local maximum or minimum occurs when the value crosses zero. There are two events, which for this example are the same. The first event (with Direction 1) corresponds to the point where the distance from the initial point is a local minimum, so that the spaceship returns to its original position. The event action is to store the time of the event in the variable tfinal and to stop the integration. The second event corresponds to a local maximum. The event action is to store the time that the spaceship is farthest from the starting position in the variable tfar. In[20]:=  |
Out[20]=
|
The first two solution components are coordinates of the body of infinitesimal mass, so plotting one against the other gives the orbit of the body. This displays one half-orbit when the spaceship is at the furthest point from the initial position. In[21]:=  |
Out[21]=
|
This displays one complete orbit when the spaceship returns to the initial position. In[22]:=  |
Out[22]=
|
Discontinuous Equations and Switching FunctionsIn many applications the function in a differential system may not be analytic or continuous everywhere. A common discontinuous problem that arises in practice involves a switching function  : In order to illustrate the difficulty in crossing a discontinuity, consider the following example [GØ84] (see also [HNW93]): Here is the input for the entire system. The switching function is assigned to the symbol event and the function defining the system depends on the sign of the switching function. In[17]:=  |
The symbol odemethod is used to indicate the numerical method that should be used for the integration. For comparison, you might want to define a different method, such as ExplicitRungeKutta, and rerun the computations in this section to see how other methods behave. This solves the system on the interval [0,1] and collects data for the mesh points of the integration using Reap and Sow. In[21]:=  |
Out[22]=
|
Here is a plot of the solution. In[23]:=  |
Out[23]=
|
Despite the fact that a solution has been obtained, it is not clear whether it has been obtained efficiently. The following example shows that the crossing of the discontinuity presents difficulties for the numerical solver. This defines a function that displays the mesh points of the integration together with the number of integration steps that are taken. In[24]:=  |
As the integration passes the discontinuity (near  in value) the integration method runs into difficulty and a large number of small steps are taken - a number of rejected steps can also sometimes be observed. In[25]:=  |
Out[25]=
|
One of the most efficient methods of crossing a discontinuity is to break the integration by restarting at the point of discontinuity. The following example shows how to use the EventLocator method to accomplish this. This numerically integrates the first part of the system up to the point of discontinuity. The switching function is given as the event. The direction of the event is restricted to a change from negative to positive. When the event is found, the solution and the time of the event are stored by the event action. In[26]:=  |
Out[28]=
|
Using the discontinuity found by the EventLocator method as a new initial condition, the integration can now be continued. This defines a system and initial condition, solves the system numerically and collects the data used for the mesh points. In[29]:=  |
Out[31]=
|
A plot of the two solutions is very similar to that obtained by solving the entire system in one go. In[32]:=  |
Out[32]=
|
Examining the mesh points it is clear that far fewer steps were taken by the method and that the problematic behavior encountered near the discontinuity has been eliminated. In[33]:=  |
Out[33]=
|
The value of the discontinuity is given as 0.6234 in [HNW93], which coincides with the value found by the EventLocator method. In this example it is possible to analytically solve the system and use a numerical method to check the value. The solution of the system up to the discontinuity can be represented in terms of Bessel functions. In[34]:=  |
Out[34]=
|
Substituting in the solution into the switching function, a local minimization confirms the value of the discontinuity. In[35]:=  |
Out[35]=
|
Avoiding Wraparound in PDEsMany evolution equations model behavior on a spatial domain which is infinite or sufficiently large to make it impractical to discretize the entire domain without using specialized discretization methods. In practice, it is often the case that it is possible to use a smaller computational domain for as long as the solution of interest remains localized. In situations where the boundaries of the computational domain are imposed by practical considerations rather than the actual model being studied, it is possible to pick boundary conditions appropriately. Using a pseudospectral method with periodic boundary conditions can make it possible to increase the extent of the computational domain because of the superb resolution of the periodic pseudospectral approximation. The drawback of periodic boundary conditions is that signals which propagate past the boundary persist on the other side of the domain, affecting the solution through wraparound. It is possible to use an absorbing layer near the boundary to minimize these effects, but it is not always possible to completely eliminate them. The sine-Gordon equation turns up in differential geometry and relativistic field theory. This example integrates the equation starting with a localized initial condition which spreads out. The periodic pseudospectral method is used for the integration. Since no absorbing layer has been instituted near the boundaries, it is most appropriate to stop the integration once wraparound becomes significant. This condition is easily detected with event location using the EventLocator method. The integration is stopped when the size of the solution at the periodic wraparound point crosses a threshold of 0.01, beyond which, the form of the wave would be affected by periodicity. In[7]:=  |
Out[7]=
|
This extracts the ending time from the InterpolatingFunction object and makes a plot of the computed solution. You can see that the integration has been stopped just as the first waves begin to reach the boundary. In[8]:=  |
Out[9]=
|
The "DiscretizedMonitorVariables" option affects the way the event is interpreted for PDEs; with the setting True, u[t,x] is replaced by a vector of discretized values. This is much more efficient because it avoids explicitly constructing the InterpolatingFunction to evaluate the event. In[10]:=  |
Out[10]=
|
Performance ComparisonThe following example constructs a table making a comparison for two different integration methods. This defines a function which returns the time it takes to compute a solution of a mildly damped pendulum equation up to the point at which the bob has momentarily been at rest 1000 times. In[36]:=  |
This uses the function to make a table comparing the different location methods for two different ODE integration methods. In[37]:=  |
Out[39]//TableForm=
|
While simple step begin/end and linear interpolation location are essentially the same low cost, the better location methods are more expensive. The default location method is particularly expensive for the explicit Runge-Kutta method because that does not yet support a continuous output formula - it therefore needs to repeatedly invoke the method with different step sizes during the local minimization. It is worth noting that, often, a significant part of the extra time for computing events arises from the need to evaluate the event functions at each time step to check for the possibility of a sign change. In[40]:=  |
Out[40]//TableForm=
|
An optimization is performed for event functions involving only the independent variable. Such events are detected automatically at initialization time. For example, this has the advantage that interpolation of the solution of the dependent variables is not carried out during at each step if the local optimization search - it is deferred until the value of the independent variable has been found. LimitationsOne limitation of the event locator method is that since the event function is only checked for sign changes over a step interval, if the event function has multiple roots in a step interval, all or some of the events may be missed. This typically only happens when the solution to the ODE varies much more slowly than the event function. When you suspect that this may have occurred, the simplest solution is to decrease the maximum step size the method can take by using the MaxStepSize option to NDSolve. More sophisticated approaches can be taken, but the best approach depends on what is being computed. An example is shown below which demonstrates the problem and shows two approaches for fixing it. This should compute the number of positive integers less than  (there are 148). However, most are missed because the method is taking large time steps because the solution x[t] is so simple. In[17]:=  |
Out[17]=
|
This restricts the maximum step size so that all the events are found. In[18]:=  |
Out[18]=
|
It is quite apparent from the nature of the example problem that if the endpoint is increased, it is likely that a smaller maximum step size may be required. Taking very small steps everywhere is quite inefficient. It is possible to introduce an adaptive time step restriction be setting up a variable which varies on the same time scale as the event function. This introduces an additional function to integrate which is the event function. With this modification and allowing the method to take as many steps is needed, it is possible to find the correct value up to t = 10 in a reasonable amount of time. In[19]:=  |
Out[19]=
|
Option summaryEventLocator optionsEventLocator method optionsEventLocationMethod optionsSettings for the EventLocationMethod option.Brent optionsOptions for event location method Brent.
| 
