Events and Discontinuities in Differential Equations

Introduction

Differential equations alone are very effective for modeling continuous behavior of systems. However, many real systems involve components that change at discrete times, possibly triggered by states of the continuous solutions. For example, in a heating system, a thermostat will switch on once the temperature reaches a certain level.

An event trigger for a differential or differential algebraic system

is a point along the solution at which a Boolean-valued event function becomes True:

When is True for , then the value of can be found numerically during integration by evaluating with an approximation of the solution and using bisection.

In general, much more robust and faster ways of numerically finding the value are available if is of the form , where is a real-valued function. Symbolic processing is automatically done on the event expression to try to get a real-valued function or simpler form if possible. For example, the Boolean event is specially processed so that an event will trigger at time , and the Boolean event is automatically processed to the real function .

For some event actions, such as discontinuity crossings, it is necessary to find a bracketing interval such that is False and is True, or for a real-valued function , the signs of and are opposite.

When an event trigger is found, typically some action is performed. In some cases, the action can affect the integration of the system, by stopping integration or even changing the values of the solution variables, so the use of events in modeling gives a great deal of flexibility for discrete actions.

In NDSolve, events augment the differential equations with WhenEvent expressions.

WhenEvent[event,action]specify an action that occurs when the event triggers it for equations in NDSolve

This section shows basic examples using WhenEvent in NDSolve. Subsequent sections in the tutorial go into more depth on topics including event detection and location, actions that modify state variables, discontinuity handling, and application examples.

These initialization commands load useful packages that have differential equations to solve and define some utility functions.
In[58]:=
Click for copyable input

A simple example is locating an event, such as when a pendulum started at a nonequilibrium 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 crosses the axis.
In[4]:=
Click for copyable input
Out[4]=

From the solution you can see that the pendulum reaches its lowest point at about . Using the 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[5]:=
Click for copyable input
Out[6]=

Event Detection and Location

Event location works in effectively two stages. In the event detection stage, for each time step taken by the numerical integration method each event is tested to see if it may have changed during the step. Once a possible event is detected, the location stage typically uses a root-finding procedure to find the actual event time .

Event Detection Methods

The simplest and fastest method for detecting whether an event occurred in a time step is to check for a change in the event function from False to True, or if it can be expressed as the root of a real function, a sign change.

The problem with just using the sign is that if a method takes large steps (such as extrapolation methods tend to) or the event function varies rapidly compared to the solution of the ODE, it is quite possible to miss events when double (or more) roots occur within a time step.

For events that can be expressed as a root of a real function, the method of using signs can be enhanced by including the time derivatives of the event function at the beginning and end of the step. For an ODE with event function , the time derivative is computed with

The values can be interpolated by a cubic Hermite polynomial that can be tested to see if it has any roots in the step interval. If necessary, the step is subdivided so that a root can be bracketed for the root-finding procedure.

Using the derivative makes the detection more robust with relatively little extra cost if does not depend on , so this is the default method for those cases. However, it is still entirely possible that events can be omitted using this method.

For some problems, it is vital to ensure that events are not missed, and NDSolve supports a method using interpolation that is very robust, but it can be significantly more expensive. This method is based on using the continuous or dense output from the method to create an interpolant of the event function with the same order as the method. The interpolant is tested for roots in the time step using interval arithmetic techniques.

The method for event detection is controlled by the option of WhenEvent.

"Sign"detect events from sign changes
"DerivativeSign"additionally use the time derivative of the event function
"Interpolation"interpolate the event function with a polynomial with the same order as the method

Settings for the option.

When the event function varies much more rapidly than the solutions of the differential equation itself, the time steps may include many crossings of an event function. This can be mitigated by including the event function in the system as a dependent variable, using the equation for the time derivative of the event shown above. This can be done automatically using the option setting "IntegrateEvent"->True in WhenEvent.

When the event is integrated, the event interpolation is even more robust, since the event function interpolant is just a component of the dense output of the time integration method. In this case, the roots of the interpolant are the event locations to the full accuracy of the method, so there is no need for using a general root-finding method. For this reason, with , the event is integrated by default.

Event Location Methods

Once an event has been detected, the next task is to refine the location of the time within the time step where the event occurs. There are several different methods that 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.

The method used to refine the location is controlled by the option of WhenEvent.

"Brent"use FindRoot with Method->"Brent" to locate the event; this is the default with the setting Automatic
"BrentBracketRoot"give the event handler a refined bracket for the root that is useful for discontinuity crossings
"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 option.

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.

Here is an example where is appropriate because the location of the event is arbitrary.

Trigger an event to stop integration when a button is clicked. If the Stop button is clicked soon enough, the integration will be terminated before .
In[7]:=
Click for copyable input
Out[7]=
In[8]:=
Click for copyable input
Out[8]=

In this example, 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, so it is most appropriate to use the or location methods. In any example where the definition of the event is heuristic or somewhat imprecise, this is the best 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 .

This computes the solution for a single period of the pendulum equation and plots the solution for that period.
In[9]:=
Click for copyable input
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]:=
Click for copyable input
Out[11]=

The linear interpolation method is sufficient for most viewing purposes, such as computing Poincaré sections. 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 , finding the location of the event using FindRoot as described in "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 and in most cases much better. You can control the accuracy and precision to which FindRoot tries to get the root of the event function using method options for the event location method. The default is to find the root to the same accuracy and precision as NDSolve is using for local error control.

When the event depends on the solutions of the differential equation, the argument needs to be approximated for within the time step. 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.

For events that are set up to handle discontinuity crossings, it is very useful for the event handler to have a bracketing interval so that it can do proper evaluations on either side of the discontinuity and ensure the crossing is done correctly. This is done with the method .

Method options can be passed to the and methods that control the solution approximation and the use of FindRoot.

option name
default value
"MaxIterations"100the maximum number of iterations to use
"AccuracyGoal"Automaticaccuracy goal setting passed to FindRoot
"PrecisionGoal"Automaticprecision goal setting passed to FindRoot
"SolutionApproximation"Automatichow to approximate the solution for evaluating the event function during the refinement process; can be Automatic or

Options for event location methods and .

If the AccuracyGoal or PrecisionGoal settings are Automatic, the value passed to FindRoot is based on the local error setting for NDSolve.

Location Method 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.
In[12]:=
Click for copyable input
This integrates the system and prints out the method used and the value of the independent variable when the integration is terminated.

Event Actions

With WhenEvent[event,action], once an event trigger point is detected and located along a solution trajectory, then the event handler evaluates the action given in the second argument of WhenEvent.

First, the solution values are approximated, typically using a dense output formula for the time integration method. The derivatives are found by evaluating the right-hand side function for ODEs or using the dense output formula for DAEs. Then (in a way similar to Block) the time-independent variable is set to and the dependent variables and their derivatives are set to the approximate values in and , and the expression is evaluated. If the value returned from the evaluation is a Rule or a string directive, then the event handler will execute the directed action.

"StopIntegration"stop integrating at ; return the solution
"RestartIntegration"restart the integration at
"CrossDiscontinuityintegrate to , extrapolate to , and restart at
"CrossSlidingDiscontinuity"integrate to , extrapolate and check sliding mode conditions, and restart at
"RemoveEvent"remove the event
y[t]->valchange the state variable y to val
d[t]->"DiscontinuitySignature"change the discontinuity signature variable d

Event action directives.

With WhenEvent[event,{action1,action2,}], is evaluated and handled completely before is evaluated, so it is possible to give a sequence of actions that will be triggered for a single event.

There are many examples of event actions with and without directives in the reference page for WhenEvent. A few examples are given here that illustrate the sequential evaluation and the handling of directives.

Collect points where the and stop integrating the first time this occurs with TemplateBox[{x}, Abs]<1.
In[14]:=
Click for copyable input
Out[14]=
Show the collected points on a plot of the solution.
In[15]:=
Click for copyable input
Out[16]=

Changing State Variables

In general, evaluation of an arbitrary action expression does not affect the solution state variables that are used by NDSolve. Using an action in WhenEvent that evaluates to a rule provides a flexible way of making state changes at a discrete event time.

A simple example of this is a model of a ball bouncing elastically on a flat surface.

Model a ball bouncing on a flat surface by reversing the velocity with a 5% loss.
In[17]:=
Click for copyable input
Out[19]=

The rule specifies that the state variable is replaced with whenever . After the change is made, the integration is restarted, since effectively a new solution trajectory is being computed.

A change for several state variables can be specified either using a single rule with a list of variables and values or multiple rules. There are subtle semantic differences between the two approaches, demonstrated by the examples below.

Swap position and velocity for an oscillator at regular time intervals.
In[20]:=
Click for copyable input
Out[22]=

In the rule , the right-hand side is evaluated and then the corresponding state variables are set. When the action is a list of rules, each rule is evaluated in turn.

Use two separate rules.
In[23]:=
Click for copyable input
Out[24]=

With the two rules, first is set to and then is set to the new value of , so the second rule actually has the effect of setting to .

For an ODE system , it is not possible to set the derivatives , since those are determined explicitly from the function. For higher-order systems that are automatically reduced to first order, this means that the highest-order derivative that appears in the equations cannot be set.

Attempt to change the derivatives.
In[25]:=
Click for copyable input
Out[27]=

When solving as a system of differential-algebraic equations , it is possible to set the derivatives. A part of the DAE solving process is to find consistent values and so that the residual is 0. NDSolve is directed to treat equations as a system of DAEs if you use the Method option setting Method->{"EquationSimplification"->"Residual"}.

Solve as DAEs.
In[28]:=
Click for copyable input
Out[30]=

When either the variables or their derivatives are changed when solving as a DAE of the form , unless the new values given satisfy , a reinitialization needs to be computed where an attempt is made to find values of the variables that were not changed, so that . Details of the procedure and examples are given in "DAE Reinitialization with Events".

When the action is not explicitly a rule, the action is evaluated after replacing all explicit instances of Rule in the action expression with a special private head that has the attribute HoldFirst, so that the states on the left-hand side do not evaluate. If the return value has this special head, then a state change is executed.

At regular intervals, make a state change with a probability and show several random realizations of the simulation.
In[31]:=
Click for copyable input
Out[33]=

Discrete variables specified with the DiscreteVariables option can only be changed through rules. An advantage of using discrete variables is that they do not increase the system size that the differential equation solver methods work withthis is particularly important for stiff methods where the Jacobian is used.

Events and Discontinuous Differential Equations

The methods built into NDSolve for solving systems of differential equations are all effectively based on an implicit assumption that the differential equation satisfies some basic continuity conditions. For ODEs, a sufficient condition is Lipschitz continuity.

When a differential equation has discontinuities, these assumptions are violated and solvers may get poor quality solutions or, in some cases, may not be able to find solutions at all. Typically what happens is that in an attempt to satisfy local error tolerances, the solvers decrease step size until the step is so small that error tolerances appear to be satisfied even with the discontinuous function values, and the solver continues beyond the discontinuity with the approximation it found.

As an example, consider the ODE

where is a parameter.

Define a function for the right-hand side of the differential equation that switches across a discontinuity.
In[34]:=
Click for copyable input
Define a function that will sample to the same values, but is hidden from symbolic evaluation.
In[37]:=
Click for copyable input
Define options specifying the method to be an explicit RungeKutta method of order 7.
In[38]:=
Click for copyable input
Solve the equation to for using the hidden definition counting the number of steps used.
In[39]:=
Click for copyable input
Out[39]=
In[40]:=
Click for copyable input
Out[40]=

The solution shown here was computed without any a priori knowledge of the discontinuity, and the reduction to the time step just prior to the discontinuity is seen. Once the discontinuity has been crossed, the time steps increase in size again. Essentially the same thing happens with the default (multistep) method, but a higher-order RungeKutta method was chosen for this example to emphasize the effect.

If the time at which a discontinuity is encountered can be identified, then an alternative way of handling the discontinuity is to step up to the discontinuity using values from the approaching side and then restart the integration using values on the other side of the discontinuity. Both the extremely small time steps and possible errors can be avoided this way.

A natural way of identifying the time at which a discontinuity is encountered is to set up an event. In the Wolfram Language, this can be done with WhenEvent. One issue with using events to detect discontinuities is that for events that depend on the solution, the solution may need to be computed past the event to correctly locate the event, but that involves the approximation of the solution across the discontinuity the event was set up to avoid in the first place. NDSolve avoids this problem by using discontinuity signature discrete variables that are only switched when the discontinuity is encountered, so that the solver in effect samples a smooth function. This is done automatically by NDSolve when it symbolically detects discontinuous functions in the differential equations, but it is possible to set up WhenEvent expressions explicitly that do the same thing; details are described following.

Solve using the explicit symbolic definition.
In[41]:=
Click for copyable input
Out[41]=

When NDSolve can see the symbolic form of the function, it is able to automatically identify the discontinuity and handle it using an event. This approach requires far fewer time steps.

This uses symbolic methods to compute a highly accurate solution and shows the error of each numerical solution at its time steps.
In[42]:=
Click for copyable input
In[47]:=
Click for copyable input
In[49]:=
Click for copyable input
Out[49]=

The error makes a jump of several orders of magnitude when the discontinuity is not handled using an event.

When the parameter is increased sufficiently, numerical methods without events are not generally able to continue past the discontinuity because when the discontinuity is first encountered, the solution on either side points back toward the discontinuity.

Attempt to solve the equation to for using the hidden definition counting the number of steps used.
In[50]:=
Click for copyable input
Out[50]=

With , numerical methods without events are not generally able to continue past the discontinuity because when the discontinuity is first encountered, the solution on either side points back toward the discontinuity.

This shows the effective vector field around the discontinuity. The discontinuity is shown in yellow, and the solution computed as far as it got is shown in red.
In[51]:=
Click for copyable input
Out[51]=
Solve using the explicit symbolic definition.
In[52]:=
Click for copyable input
Out[52]=
Plot the solution, showing the time steps as points and the discontinuity in yellow.
In[53]:=
Click for copyable input
Out[53]=

With the symbolic detection of the discontinuity, even though the Wolfram Language conditions for unique existence of the solution are violated, NDSolve is able to compute a continuation known as Filippov [F88] sliding mode in this situation. When the vector field points toward the discontinuity on both sides, the continuation effectively slides along the discontinuity surface until the field on one side or the other points away. This can be seen in the plot above where the solution is directly on the circle for some time, but eventually deviates.

The following sections of this tutorial explain the Filippov continuation and how to use discontinuity signatures.

Sliding Mode

When a discontinuity surface is encountered and the vector field on either side of the surface points back toward the solution, some sort of continuation method is necessary. Without a continuation, a numerical method will typically fail as it tries to take steps that cross and then recross the discontinuity; in some cases the numerical solution will display artificial chattering.

An appropriate continuation defined by Filippov [F88] effectively takes a convex combination of the components of the vector fields pointing along the discontinuity surface (or orthogonal to the gradient vector of a function for which the surface is a level curve).

114.gif

Entering sliding mode.

In the preceding figure, the discontinuity surface is defined by and the vector fields for and are and , respectively. Then on the surface use , where is chosen so that is normal to the gradient , in particular . The solution can enter sliding mode either from above or below, as shown in the left and right parts of the figure, respectively.

The preceding definition is given for an autonomous system with no explicit dependence on the independent variable . The same ideas apply for a nonautonomous system by adding a variable and considering the autonomous system .

A simple example is provided by the vector field that has a discontinuity surface . For this vector field , , and ; therefore , so and the sliding motion is governed by .

Plot the vector field.
In[54]:=
Click for copyable input
Out[54]=
This attempts to solve the problem without any handling of the discontinuity surface.
In[57]:=
Click for copyable input
Out[57]=

The default method gets stuck because it eventually gets to tiny steps that chatter along the discontinuity.

This attempts to solve the problem using an extrapolation method.
In[58]:=
Click for copyable input
In[60]:=
Click for copyable input
Out[60]=

The solution shown above is not correct; the larger oscillation or chattering of the solution near the sliding surface is an artifact of the numerical method.

Solve the problem with automatic discontinuity handling.
In[61]:=
Click for copyable input
In[62]:=
Click for copyable input
Out[62]=

When the discontinuous function is explicit in the equations, NDSolve can handle the discontinuity automatically by using the Filippov continuation.

The Filippov continuation is a natural choice since solving a sequence of systems with smooth functions approaching the discontinuous one approaches the sliding mode solution.

Define a smooth function that approaches as the parameter goes to zero.
In[63]:=
Click for copyable input
Out[64]=
Manipulate the solution over a range of .

In general, a trajectory will be trapped in sliding mode until one of the vector fields starts to point away. When this happens, the solution can be continued above or below the discontinuity.

145.gif

Exiting sliding mode.

In the figure above, the discontinuity surface is defined by and the vector fields for and are and , respectively. Then on the surface use , where is chosen so that is normal to the gradient , in particular . Beyond the point shown in both sides of the figure, and both have the same sign. On the left because and , and on the right because and , and in either case the trajectory exits sliding mode. These points can be found by numerically finding either the root or as appropriate.

Consider the vector field . Once , the vector field shown following starts to point down again, and the trajectory can escape sliding mode.

Solve the system using sliding mode and show the solution with the vector field.
In[66]:=
Click for copyable input
In[67]:=
Click for copyable input
Out[67]=

It may occur that a solution enters and exits sliding mode repeatedly.

Consider the vector field . For , this is just an oscillator with negative damping. For , all trajectories approach the discontinuity surface . On , if is not too large, there is a set of points where there is sliding motion. The sliding motion in effect stabilizes the unbounded growth that would occur with the negative damping, and so there is a periodic orbit.

This sets up a dynamic plot that allows you to move the initial condition.

All of the examples shown so far are in two dimensions, since that is easier to visualize. However, the continuation works in any number of dimensions, and solutions can arise where the solution slides on more than one discontinuity surface at once.

Consider the two discontinuity surfaces in 3D defined by the surface of a sphere and the plane and the vector field given by .

This solves the system with two discontinuities.
In[69]:=
Click for copyable input
Out[69]=
Out[70]=
This shows the trajectory in 3D.
In[71]:=
Click for copyable input
In[73]:=
Click for copyable input
Out[74]=

The solution first slides along the surface of the sphere, then along both surfaces, and eventually just along the plane.

Discontinuity Signature

For functions with complicated definitions such that symbolic processing of discontinuities is not possible, it is possible to set up manually the same sort of discontinuity handling that NDSolve does automatically for discontinuous functions.

Discontinuity signature variables are set up by including a WhenEvent expression where the event is on the discontinuity surface, and the action is rule setting the variable to and using the DiscreteVariables option to specify the range of the variable to be either or .

WhenEvent[e0,s[t]->"DiscontinuitySignature]specify that the discrete variable holds the discontinuity signature of the discontinuity surface defined by e such that is effectively Sign[e]
DiscreteVariables->{s[t]{-1,0,1}}indicate that is a discrete variable with possible values , , and ; sliding mode solutions will be considered
DiscreteVariables->{s[t]{-1,1}}indicate that is a discrete variable with possible values and ; sliding mode solutions will not be considered

Setting up a discontinuity signature variable .

Consider the first example in the previous section on sliding mode with the vector field .

Solve the system using an equivalent formulation with a discontinuity signature variable.
In[75]:=
Click for copyable input
Out[75]=
Plot the solution and the values of the discontinuity signature variable .
In[76]:=
Click for copyable input
Out[76]=

The variable switches from to 0 when the discontinuity is encountered and the solution goes into sliding mode. Sliding mode is indicated by a discontinuity signature of 0, since ideally the solution is exactly on the discontinuity surface and so is 0. In practice, the actual solution may not satisfy this due to numerical error. When in sliding mode, NDSolve tries to keep as close to the surface as possible.

If you know that a discontinuity will not lead to sliding mode, the needed computations can be done less expensively if you exclude 0 from the range of the discontinuity signature variable. However, if you do this and there could be sliding mode on the discontinuity, the solution will fail or be incorrect.

Attempt to solve the system using an equivalent formulation with a discontinuity signature variable.
In[77]:=
Click for copyable input
Out[77]=
The solution is incorrect.
In[78]:=
Click for copyable input
Out[78]=

An autonomous system defined by

can be represented by , where is set by the Wolfram Language expression WhenEvent[e[x[t]]>0,s[t]->"DiscontinuitySignature"] in NDSolve.

The time derivative of the event function, is important in determining what happens when at some time . There are six cases where will be changed when with WhenEvent[e[x[t]]>0,s[t]->"DiscontinuitySignature"].

206.gif

Crossing from below and above.

  • Case 1: (Cross from below) and for with at . The solution trajectory simply crosses to and is set to .
  • Case 2: (Cross from above) and for with at .The solution trajectory simply crosses to and is set to .

223.gif

Entering sliding mode from above and below.

  • Case 3: (Enter sliding mode from above) and for with at . Since is decreasing to zero it also holds that at . The derivative is negative above and positive below the discontinuity, so the trajectory is trapped on the discontinuity. This is referred to as sliding mode and the value of is set to 0. The effective vector field is computed using the Filippov continuation , where . Note that since and .
  • Case 4: (Enter sliding mode from below) and for with at . This is essentially the same as case 3 and the solution enters sliding mode with set to 0.

246.gif

Exiting sliding mode to above and below.

  • Case 5: (Exit sliding mode to above) , , , and for and at with . The solution is no longer in sliding mode and is set to 1. Note that the event that triggers this change is effectively a zero crossing of .
  • Case 6: (Exit sliding mode to below) , , , and for and at with . The solution is no longer in sliding mode and is set to . Note that the event that triggers this change is effectively a zero crossing of .

To see all of these cases in an actual example, consider the system of ODEs

where and are constants that can be defined to see the different cases.

Define a function that gives the differential equations given definitions fa for and fb for .
In[79]:=
Click for copyable input
Define a function solver that, given function definitions and an initial condition for , solves the system using NDSolve with as a discontinuity signature variable.
In[80]:=
Click for copyable input
Define a function that plots the solution superimposed on the vector field and a plot of the discontinuity signature variable along with if there is sliding mode.
In[81]:=
Click for copyable input
Case 1: Cross from below with and .
In[82]:=
Click for copyable input
Out[82]=
Case 2: Cross from above with and .
In[83]:=
Click for copyable input
Out[83]=
Case 3: Enter sliding mode from above with and .
In[84]:=
Click for copyable input
Out[84]=
Case 4: Enter sliding mode from below with and .
In[85]:=
Click for copyable input
Out[85]=
Case 5: Exit sliding mode to above and .
In[86]:=
Click for copyable input
Out[86]=
Case 6: Exit sliding mode to below with and .
In[87]:=
Click for copyable input
Out[87]=

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.
In[88]:=
Click for copyable input
Out[88]=
This plots the solution and highlights the initial and final points (green and red) by encircling them.
In[89]:=
Click for copyable input
Out[89]=

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.
In[90]:=
Click for copyable input
Out[90]=

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 .
In[91]:=
Click for copyable input
This uses the function to compute the period at .
In[92]:=
Click for copyable input
Out[92]=
This shows a log-log plot of the period as a function of .
In[93]:=
Click for copyable input
Out[93]=

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énonHeiles System

Define the HénonHeiles system that models stellar motion in a galaxy.

This gets the HénonHeiles system from the package.
In[94]:=
Click for copyable input
Out[95]=

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.

Specify a WhenEvent that defines the Poincaré section. The event action is to use Sow on the values of and .
In[96]:=
Click for copyable input
This integrates the HénonHeiles system using a fourth-order explicit RungeKutta method with fixed step size of 0.25. The event action is to use Sow on the values of and .
In[97]:=
Click for copyable input
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[98]:=
Click for copyable input
Out[99]=

Since the HénonHeiles system is Hamiltonian, a symplectic method gives much better qualitative results for this example.

This integrates the HénonHeiles system using a fourth-order symplectic partitioned RungeKutta method with fixed step size of 0.25. The event action is to use Sow on the values of and .
In[100]:=
Click for copyable input
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[101]:=
Click for copyable input
Out[102]=

The ABC Flow

This loads an example problem of the ArnoldBeltramiChildress (ABC) flow that is used to model chaos in laminar flows of the three-dimensional Euler equations.
In[103]:=
Click for copyable input
This defines a splitting , of the system by setting some of the right-hand side components to zero.
In[107]:=
Click for copyable input
Out[107]=
In[108]:=
Click for copyable input
Out[108]=
This defines the implicit midpoint method.
In[109]:=
Click for copyable input
This constructs a second-order splitting method that retains volume and reversing symmetries.
In[110]:=
Click for copyable input
This defines a function that gives the Poincaré section for a particular initial condition.
In[111]:=
Click for copyable input
This finds the Poincaré sections for several different initial conditions and flattens them together into a single list of points.
In[112]:=
Click for copyable input
Out[113]=

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 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.
In[114]:=
Click for copyable input
This defines the function that runs the bouncing ball simulation for a given reflection ratio, ramp, and starting position.
In[115]:=
Click for copyable input
This is the example that is done in [SGT03].
In[116]:=
Click for copyable input
Out[117]=
The ramp is now defined to be a quarter circle.
In[118]:=
Click for copyable input
Out[119]=
This adds a slight waviness to the ramp.
In[120]:=
Click for copyable input
Out[121]=

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 as described in "Pseudospectral Derivatives" 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 WhenEvent.

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[122]:=
Click for copyable input
Out[123]=
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[124]:=
Click for copyable input
Out[125]=
The option affects the way the event is interpreted for PDEs; with the setting True, 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[126]:=
Click for copyable input
Out[127]=

Three-Body Problem

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.
In[128]:=
Click for copyable input
The event function is the derivative of the distance from the initial conditions. A local maximum or minimum occurs when the value crosses zero.
In[133]:=
Click for copyable input
There are two events. The first event, where the distance derivative becomes negative, corresponds to a local maximum in distance. The second event corresponds to the point where the distance from the initial point is a local minimum, so that the spaceship returns to its original position. Note that Evaluate is necessary in the first argument of WhenEvent so that the variable substitutions get done correctly.
In[134]:=
Click for copyable input
Out[134]=

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 complete orbit when the spaceship returns to the initial position. The far point is highlighted with a point.
In[135]:=
Click for copyable input
Out[135]=

Infection Model with Delay

The following system models an infectious disease (see [HNW93], [ST00], and [ST01]).
In[136]:=
Click for copyable input
Define a function that makes a WhenEvent expression to capture the local maximum of the ^(th) component. A separate tag for Sow is used to distinguish the components.
In[138]:=
Click for copyable input
Solve the system with a WhenEvent for the maximum of each component.
In[139]:=
Click for copyable input
Out[139]=
Display the local maxima together with the solution components.
In[140]:=
Click for copyable input
Out[143]=

Discontinuity Signature and C Code

For symbolically specified differential equations, NDSolve will automatically scan the equations for discontinuous functions and set up events for handling the discontinuities appropriately. If a differential system is defined by evaluating a function that it is not possible to symbolically scan, then it may be necessary to set up discontinuity signature variables and events to handle discontinuities that may be embedded in the definition.

A function defined in C code that is compiled to link to the Wolfram Language with LibraryFunction is a good example of a function that the Wolfram Language cannot symbolically scan.

C source code string for the right-hand side of a differential system.
In[144]:=
Click for copyable input
Use CreateLibrary from the CCompilerDriver package to make a library.
In[145]:=
Click for copyable input
In[146]:=
Click for copyable input
Out[146]=
Load the library into a LibraryFunction.
In[147]:=
Click for copyable input
Out[147]=
Solve the system for with the switching on the circle .
In[148]:=
Click for copyable input
Out[148]=
Show a plot of the solution with blue for , green for , and red for sliding mode.
In[149]:=
Click for copyable input
Out[149]=
With automatic discontinuity handling, an equivalent solution can be found with symbolic input.
In[150]:=
Click for copyable input
Out[151]=