"EventLocator" Method for NDSolve

Introduction	Limitations
Event Location Methods	Option Summary
Examples

Introduction

It 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 that 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:

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 y[t] crosses the axis:

From the solution you can see that the pendulum reaches its lowest point at about . 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:

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 earlier. 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

is detected for a damped pendulum:

Note that in the 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 that are useful for collecting data when you do not, 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 you to collect data for each instance of an event:

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 that use FindRoot to control tolerances for finding the root.

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 that, when clicked, changes the value of the variable stop to True from its initialized value of False:

This integrates the pendulum equation up until the button is clicked (or the system runs out of memory):

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 button is clicked. The number of complete swings of the pendulum is kept track of during the integration:

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 are based on the values at the initial condition and need to be consistent throughout the integration.

The "EventCondition" option of "EventLocator" allows you to specify additional Boolean conditions that need to be satisfied for an event to be tested. It is advantageous to use this instead of a Boolean event when possible because the root finding process can be done more efficiently.

This stops the integration of a damped pendulum at the first time that

once the decay has reduced the energy integral to

This makes a plot of the solution (black), the derivative (blue), and the energy integral (green). The energy threshold is shown in red:

The Method option of "EventLocator" allows the specification of the numerical method to use in the integration.

Event Location Methods

The "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 to 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" Methods

The crudest methods are appropriate for when the exact position of the event location does not 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 you can time the click 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" Method

When 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.

Denote the event function values at successive mesh points of the numerical integration:

Linear interpolation gives

A linear approximation of the event time is then:

Linear interpolation could also be used to approximate the solution at the event time. However, since derivative values and are available at the mesh points, a better approximation of the solution at the event can be computed cheaply using cubic Hermite interpolation as

for suitably defined interpolation weights:

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:

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:

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 Method

The 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 that 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 that 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 that 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.

Comparison

This 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:

Examples

Falling Body

This 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:

This plots the solution and highlights the initial and final points (green and red) by encircling them:

Period of the Van der Pol Oscillator

The 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

up to the point where the variable

reaches its initial value and direction:

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 that returns the period as a function of .

This defines a function that returns the period as a function of

This uses the function to compute the period at

Of course, it is easy to generalize this method to any system with periodic solutions.

Poincaré Sections

Using Poincaré sections is a useful technique for visualizing the solutions of high-dimensional differential systems.

For an interactive graphical interface see the package EquationTrekker.

The Hénon–Heiles System

Define the Hénon–Heiles system that models stellar motion in a galaxy.

This gets the Hénon–Heiles system from the NDSolveProblems package:

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 InterpolatingFunction by specifying the output variables list (in the second argument to NDSolve) as empty, or {}. This means that NDSolve will produce no InterpolatingFunction 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 to 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:

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 use Sow on the values of

and

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:

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 use Sow on the values of

and

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:

The ABC Flow

This loads an example problem of the Arnold–Beltrami–Childress (ABC) flow that is used to model chaos in laminar flows of the three-dimensional Euler equations:

This defines a splitting Y1, Y2 of the system by setting some of the right-hand side components to zero:

This defines the implicit midpoint method:

This constructs a second-order splitting method that retains volume and reversing symmetries:

This defines a function that gives the Poincaré section for a particular initial condition:

This finds the Poincaré sections for several different initial conditions and flattens them together into a single list of points:

Bouncing Ball

This example is a generalization of an example in [SGT03]. It models a ball bouncing down a ramp with a given profile. The example is good for demonstrating how you can use multiple invocations of NDSolve with event location to model some behavior.

This defines a function that computes the solution from one bounce to the next. The solution is computed until the next time the path intersects the ramp:

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

of energy is left after a bounce:

This defines the function that runs the bouncing ball simulation for a given reflection ratio, ramp, and starting position:

This is the example that is done in [SGT03]:

The ramp is now defined to be a quarter circle:

This adds a slight waviness to the ramp:

Event Direction

Ordinary Differential Equation

This example illustrates the solution of the restricted three-body problem, a standard nonstiff test system of four equations. The example traces the path of a spaceship traveling 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 traveling 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:

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:

This displays one complete orbit when the spaceship returns to the initial position:

Delay Differential Equation

The following system models an infectious disease (see [HNW93], [ST00], and [ST01]):

Collect the data for a local maximum of each component as the integration proceeds. A separate tag for Sow and Reap is used to distinguish the components:

Display the local maxima together with the solution components:

Discontinuous Equations and Switching Functions

In 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:

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:

Here is a plot of the solution:

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:

As the integration passes the discontinuity (near 0.6 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:

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:

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:

A plot of the two solutions is very similar to that obtained by solving the entire system at once:

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:

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 and gamma functions:

Substituting in the solution into the switching function, a local minimization confirms the value of the discontinuity:

Avoiding Wraparound in PDEs

Many evolution equations model behavior on a spatial domain that 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 that 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 that 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:

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:

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:

Performance Comparison

The following example constructs a table making a comparison for two different integration methods.

This defines a function that 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:

This uses the function to make a table comparing the different location methods for two different ODE integration methods:

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 it 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.

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 at each step of the local optimization search—it is deferred until the value of the independent variable has been found.

Limitations

One 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 follows that 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:

This restricts the maximum step size so that all the events are found:

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 by setting up a variable that varies on the same time scale as the event function.

This introduces an additional function to integrate: the event function. With this modification and allowing the method to take as many steps as needed, it is possible to find the correct value up to

in a reasonable amount of time:

Option Summary

"EventLocator" Options

option name	default value
"Direction"	All	the direction of zero crossing to allow for the event; 1 means from negative to positive, -1 means from positive to negative, and All includes both directions
"Event"	None	an expression that defines the event; an event occurs at points where substituting the numerical values of the problem variables makes the expression equal to zero
"EventAction"	Throw[Null, "StopIntegration"]	what to do when an event occurs: problem variables are substituted with their numerical values at the event; in general, you need to use RuleDelayed () to prevent the option from being evaluated except with numerical values
"EventLocationMethod"	Automatic	the method to use for refining the location of a given event
"Method"	Automatic	the method to use for integrating the system of ODEs

"EventLocator" method options.

"EventLocationMethod" Options

"Brent"	use FindRoot with Method->"Brent" to locate the event; this is the default with the setting Automatic
"LinearInterpolation"	locate the event time using linear interpolation; cubic Hermite interpolation is then used to find the solution at the event time
"StepBegin"	the event is given by the solution at the beginning of the step
"StepEnd"	the event is given by the solution at the end of the step

Settings for the "EventLocationMethod" option.

"Brent" Options

option name	default value
"MaxIterations"	100	the maximum number of iterations to use for locating an event within a step of the method
"AccuracyGoal"	Automatic	accuracy goal setting passed to FindRoot; if Automatic, the value passed to FindRoot is based on the local error setting for NDSolve
"PrecisionGoal"	Automatic	precision goal setting passed to FindRoot; if Automatic, the value passed to FindRoot is based on the local error setting for NDSolve
"SolutionApproximation"	Automatic	how to approximate the solution for evaluating the event function during the refinement process; can be Automatic or "CubicHermiteInterpolation"

Options for event location method "Brent".

Top