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 (Subscript[t, a]) along the solution at which a Boolean-valued event function becomes (RefLink[True, paclet:ref/True]):

When b is (RefLink[True, paclet:ref/True]) for t≥(Subscript[t, a]), then the value of (Subscript[t, a]) can be found numerically during integration by evaluating b with an approximation of the solution and using bisection.

In general, much more robust and faster ways of numerically finding the value (Subscript[t, a]) are available if b is of the form g(t,X(t),X'(t))==0, where g 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 b==(t==(Subscript[t, 0])) is specially processed so that an event will trigger at time (Subscript[t, a])==(Subscript[t, 0]), and the Boolean event b==(sin(t)==1/2) is automatically processed to the real function g==sin(t)-1/2 .

For some event actions, such as discontinuity crossings, it is necessary to find a bracketing interval (Subscript[t, a-])<(Subscript[t, a])<(Subscript[t, a+]) such that b((Subscript[t, a-]),X((Subscript[t, a-])),X'((Subscript[t, a-]))) is (RefLink[False, paclet:ref/False]) and b((Subscript[t, a+]),X((Subscript[t, a+])),X'((Subscript[t, a+]))) is (RefLink[True, paclet:ref/True]), or for a real-valued function g, the signs of g((Subscript[t, a-]),X((Subscript[t, a-])),X'((Subscript[t, a-]))) and g((Subscript[t, a+]),X((Subscript[t, a+])),X'((Subscript[t, a+]))) 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 (RefLink[NDSolve, paclet:ref/NDSolve]), events augment the differential equations with (RefLink[WhenEvent, paclet:ref/WhenEvent]) expressions.

This section shows basic examples using (RefLink[WhenEvent, paclet:ref/WhenEvent]) in (RefLink[NDSolve, paclet:ref/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:

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];

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

sol = NDSolve[{y''[t] + Sin[y[t]] == 0, y'[0] == 0, y[0] == 1, WhenEvent[y[t]==0,"StopIntegration"]},y,{t,0,10}]

From the solution you can see that the pendulum reaches its lowest point y(t)==0 at about t==1.675. Using the InterpolatingFunctionAnatomy package, it is possible to extract the value from the (RefLink[InterpolatingFunction, paclet:ref/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:

end = InterpolatingFunctionDomain[First[y /. sol]][[1,-1]];
Plot[Evaluate[{y[t], y'[t]} /. First[sol]], {t,0,end}, PlotStyle->{{Black},{Blue}}]

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 (Subscript[t, a]).

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 (RefLink[False, paclet:ref/False]) to (RefLink[True, paclet:ref/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 X'(t)==f(t,X(t)) with event function g(t,X(t)), 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 g does not depend on X', 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 (RefLink[NDSolve, paclet:ref/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 "EventDetection" option of (RefLink[WhenEvent, paclet:ref/WhenEvent]).

Settings for the "DetectionMethod" 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"->(RefLink[True, paclet:ref/True]) in (RefLink[WhenEvent, paclet:ref/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 "DetectionMethod"->"Interpolation", 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 "LocationMethod" option of (RefLink[WhenEvent, paclet:ref/WhenEvent]).

Settings for the "LocationMethod" 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 "StepEnd" 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 t==100:

Button["Stop",stop = True]

stop = False;usol = NDSolveValue[{D[u[t,x],t,t]==D[u[t,x],x,x], u[0,x]==Exp[-10 x^2],Derivative[1,0][u][0,x]==0,u[t,-10]==u[t,10],WhenEvent[stop,"StopIntegration"]},u,{t,0,100},{x,-10,10}, Method->"StiffnessSwitching"]

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 "StepBegin" or "StepEnd" 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 (Subscript[f, n])==f((Subscript[t, n]),(Subscript[y, n]))  and (Subscript[f, n+1])==f((Subscript[t, n+1]),(Subscript[y, n+1])) 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:

lsol = NDSolveValue[{y''[t] + Sin[y[t]] == 0, y[0] == 3, y'[0] == 0,WhenEvent[y'[t]<0,{lend=t,"StopIntegration"}, "LocationMethod"->"LinearInterpolation"]},y,{t,0,∞}, Method->"ExplicitRungeKutta"];
Plot[{lsol[t],lsol'[t]},{t,0,lend}]

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:

Plot[lsol'[t], {t,end*(1 - .001),end*(1+.001)}, PlotStyle->Blue, PlotRange->All,Epilog->{PointSize[.025],Point[{end,0}]}]

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 "Brent", finding the location of the event using (RefLink[FindRoot, paclet:ref/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 (RefLink[FindRoot, paclet:ref/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 (RefLink[NDSolve, paclet:ref/NDSolve]) is using for local error control.

When the event depends on the solutions x(t) of the differential equation, the argument x(t) needs to be approximated for t 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 (RefLink[FindRoot, paclet:ref/FindRoot]) on the (RefLink[InterpolatingFunction, paclet:ref/InterpolatingFunction]) object returned from (RefLink[NDSolve, paclet:ref/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 (Subscript[t, a-])<(Subscript[t, a])<(Subscript[t, a+]) 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 "BrentBracketRoot".

Method options can be passed to the "Brent" and "BrentBracketRoot" methods that control the solution approximation and the use of (RefLink[FindRoot, paclet:ref/FindRoot]).

Options for event location methods "Brent" and "BrentBracketRoot".

If the (RefLink[AccuracyGoal, paclet:ref/AccuracyGoal]) or (RefLink[PrecisionGoal, paclet:ref/PrecisionGoal]) settings are (RefLink[Automatic, paclet:ref/Automatic]), the value passed to (RefLink[FindRoot, paclet:ref/FindRoot]) is based on the local error setting for (RefLink[NDSolve, paclet:ref/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:

eventmethods = {"StepBegin","StepEnd","LinearInterpolation", Automatic};

This integrates the system and prints out the method used and the value of the independent variable when the integration is terminated:

Map[
NDSolve[{y''[t] + Sin[y[t]] == 0, y[0] == 3, y'[0] == 0,WhenEvent[y'[t]<0,{Print[#,": t = ",t," y'[t] = ", y'[t]],"StopIntegration"},"LocationMethod"->#]},y,{t,0,∞}, Method->"ExplicitRungeKutta"]&,
eventmethods
];

Event Actions

With (RefLink[WhenEvent, paclet:ref/WhenEvent])[event,action], once an event trigger point (Subscript[t, a]) is detected and located along a solution trajectory, then the event handler evaluates the action given in the second argument of (RefLink[WhenEvent, paclet:ref/WhenEvent]).

First, the solution values X((Subscript[t, a])) are approximated, typically using a dense output formula for the time integration method. The derivatives X'((Subscript[t, a])) 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 (RefLink[Block, paclet:ref/Block])) the time-independent variable is set to (Subscript[t, a]) and the dependent variables and their derivatives are set to the approximate values in X((Subscript[t, a])) and X'((Subscript[t, a])), and the expression is evaluated. If the value returned from the evaluation is a (RefLink[Rule, paclet:ref/Rule]) or a string directive, then the event handler will execute the directed action.

Event action directives.

With (RefLink[WhenEvent, paclet:ref/WhenEvent])[event,{(Subscript[action, 1]),(Subscript[action, 2]),…}], (Subscript[action, i]) is evaluated and handled completely before (Subscript[action, i+1]) 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 (RefLink[WhenEvent, paclet:ref/WhenEvent]). A few examples are given here that illustrate the sequential evaluation and the handling of directives.

Collect points where the x==y and stop integrating the first time this occurs with (Abs[x])<1:

{sol, data} = Reap[NDSolve[{x'[t]==y[t],y'[t]==-t x[t],x[0]==1,y[0]==1,WhenEvent[x[t]==y[t],Sow[{t,x[t]}];If[Abs[x[t]]<1,"StopIntegration"]]},{x,y},{t,∞}]]

Show the collected points on a plot of the solution:

tend= data[[1,-1,1]];
Plot[Evaluate[{x[t],y[t]} /.sol],{t,0,tend},Epilog->{PointSize[0.05],Red,Apply[Point,data]}]

Changing State Variables

In general, evaluation of an arbitrary action expression does not affect the solution state variables that are used by (RefLink[NDSolve, paclet:ref/NDSolve]). Using an action in (RefLink[WhenEvent, paclet:ref/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:

g = 9.8; r = .95;
bb = NDSolveValue[{x''[t]==-g,x[0]==1,x'[0]==0,WhenEvent[x[t]==0,x'[t]->-r x'[t]]},x,{t,0,5}];
Plot[bb[t],{t,0,5}]

The rule x'[t]->-r x[t] specifies that the state variable x'[t] is replaced with -r x'[t] whenever x[t]==0. 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:

Δt=3;T = 20;
sol = First[NDSolve[{x'[t]==y[t],y'[t]+  y[t]+2 x[t]==0,x[0]==0,y[0]==1,WhenEvent[Mod[t,Δt]==0,{x[t],y[t]}->{4 y[t],x[t]}]},{x,y},{t,0,T}]];
ParametricPlot[Evaluate[{x[t],y[t]} /. sol],{t,0,T},PlotRange->All]

In the rule {x[t],y[t]}->{y[t],x[t]}, 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:

sol2 = First[NDSolve[{x'[t]==y[t],y'[t]+  y[t]+2 x[t]==0,x[0]==0,y[0]==1,WhenEvent[Mod[t,Δt]==0,{x[t]->4 y[t],y[t]->x[t]}]},{x,y},{t,0,T}]];
ParametricPlot[Evaluate[{x[t],y[t]} /. sol2],{t,0,T},PlotRange->All]

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

For an ODE system x'(t)==f(t,x(t)), it is not possible to set the derivatives x'[t], 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:

Δt=3;T = 20;
dsol = First[NDSolve[{x'[t]==y[t],y'[t]+  y[t]+2 x[t]==0,x[0]==0,y[0]==1,WhenEvent[Mod[t,Δt]==0,{x'[t],y'[t]}->{4 y[t],x[t]}]},{x,y},{t,0,T}]];
ParametricPlot[Evaluate[{x[t],y[t]} /. dsol],{t,0,T},PlotRange->All]

When solving as a system of differential-algebraic equations g(t,x(t),x'(t))==0, it is possible to set the derivatives. A part of the DAE solving process is to find consistent values x(t) and x'(t) so that the residual g is 0. (RefLink[NDSolve, paclet:ref/NDSolve]) is directed to treat equations as a system of DAEs if you use the (RefLink[Method, paclet:ref/Method]) option setting (RefLink[Method, paclet:ref/Method])->{"EquationSimplification"->"Residual"}.

Solve as DAEs:

Δt=3;T = 20;
daesol = First[NDSolve[{x'[t]==y[t],y'[t]+  y[t]+2 x[t]==0,x[0]==0,y[0]==1,WhenEvent[Mod[t,Δt]==0,{x'[t],y'[t]}->{4 y[t],x[t]}]},{x,y},{t,0,T}, Method->{"EquationSimplification"->"Residual"}]];
ParametricPlot[Evaluate[{x[t],y[t]} /. daesol],{t,0,T},PlotRange->All]

When either the variables or their derivatives are changed when solving as a DAE of the form g(t,x(t),x'(t))==0, unless the new values given satisfy g==0, a reinitialization needs to be computed where an attempt is made to find values of the variables that were not changed, so that g==0. 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 (RefLink[Rule, paclet:ref/Rule]) in the action expression with a special private head that has the attribute (RefLink[HoldFirst, paclet:ref/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 (1/3) probability and show several random realizations of the simulation:

Δt=1;
xsols = Table[NDSolveValue[{x''[t]+x'[t] + x[t]==Cos[t],x[0]==x'[0]==0,WhenEvent[Mod[t,Δt]==0,If[RandomChoice[{1/3,2/3}->{True,False}],x[t]->-x[t]]]},x[t],{t,0,10}],{10}];
Plot[Evaluate[xsols],{t,0,10}]

Discrete variables specified with the (RefLink[DiscreteVariables, paclet:ref/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 with—this is particularly important for stiff methods where the Jacobian is used.

Events and Discontinuous Differential Equations

The methods built into (RefLink[NDSolve, paclet:ref/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 a is a parameter.

Define a function f for the right-hand side of the differential equation that switches across a discontinuity:

f1[a_][t_,x_] = (t^2) + 2 (x^2);
f2[a_][t_,x_]=2(t^2) + 3 (x^2)-a;
f[a_][t_,x_] = If[(((x+(3/20)))^2)+(((t+(1/20)))^2)≤1,f1[a][t,x],f2[a][t,x]];

Define a function that will sample to the same values, but is hidden from symbolic evaluation:

fb[a_][t_?NumberQ,x_]:= f[a][t,x];

Define options specifying the method to be an explicit Runge–Kutta method of order 7:

opts = {Method->{"TimeIntegration"->{"ExplicitRungeKutta", "DifferenceOrder"->7}}};

Solve the equation to t==1 for a==2 using the hidden definition counting the number of steps used:

Block[{s =0},bsol = NDSolveValue[{x'[t]==fb[2][t,x[t]],x[0]==3/10},x,{t,1}, StepMonitor:>s++,opts];s]

ListLinePlot[bsol, Mesh->All,PlotRange->All]

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 Runge–Kutta 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 (RefLink[WhenEvent, paclet:ref/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. (RefLink[NDSolve, paclet:ref/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 (RefLink[NDSolve, paclet:ref/NDSolve]) when it symbolically detects discontinuous functions in the differential equations, but it is possible to set up (RefLink[WhenEvent, paclet:ref/WhenEvent]) expressions explicitly that do the same thing; details are described following.

Solve using the explicit symbolic definition:

Block[{s =0},dsol = NDSolveValue[{x'[t]==f[2][t,x[t]],x[0]==3/10},x,{t,1}, StepMonitor:>s++, opts];s]

When (RefLink[NDSolve, paclet:ref/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:

x1[t_] = First[x[t] /. DSolve[{x'[t]==f1[2][t,x[t]], x[0] == 3/10}, x, t]];
switchpoint = t /. FindRoot[(((x+(3/20)))^2)+(((t+(1/20)))^2)==1 /. {x->x1[t]}, {t, .6}, WorkingPrecision->100];
x2[t_] = Block[{a = 2}, First[x[t] /. DSolve[{x'[t] ==f2[2][t,x[t]], x[s] == x1[s]}, x, t] /. s->switchpoint]];
SetAttributes[x100, Listable];
x100[t_] := Re[Piecewise[{{x1[t], t≤ switchpoint}}, x2[t]]];

bsteps = Part[InterpolatingFunctionCoordinates[bsol],1,2;;-1];
dsteps = Part[InterpolatingFunctionCoordinates[dsol],1,2;;-1];

ListLogPlot[{Map[{#, Abs[bsol[#] - x100[#]]}&, bsteps], Map[{#, Abs[dsol[#] - x100[#]]}&, dsteps]}, PlotRange->All,Joined->True]

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

When the parameter a 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 t=1 for a=5 using the hidden definition counting the number of steps used:

Block[{s =0},bsol = NDSolveValue[{x'[t]==fb[5][t,x[t]],x[0]==3/10},x,{t,1}, StepMonitor:>s++,opts];s]

With a=5, 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:

Show[VectorPlot[{1,f[5][t,x]},{t,.5,.7},{x,.5,.7},Prolog->{Thickness[.015], Yellow, Circle[{-1/20, -3/20}]}, FrameLabel->{t,x}], Plot[bsol[t],{t,.5,0.6235830505886174}, PlotStyle->Red]]

Solve using the explicit symbolic definition:

Block[{s =0},dsol = NDSolveValue[{x'[t]==f[5][t,x[t]],x[0]==3/10},x,{t,1}, StepMonitor:>s++,opts];s]

Plot the solution, showing the time steps as points and the discontinuity in yellow:

ListLinePlot[dsol, Mesh->All,PlotRange->All,Prolog->{Thickness[.015], Yellow, Circle[{-1/20, -3/20}]},AxesLabel->{t,x}]

With the symbolic detection of the discontinuity, even though the Wolfram Language conditions for unique existence of the solution are violated, (RefLink[NDSolve, paclet:ref/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).

Entering sliding mode.

In the preceding figure, the discontinuity surface is defined by e==0 and the vector fields for e>0 and e<0 are (Subscript[f, 1]) and (Subscript[f, -1]), respectively. Then on the surface ϕ use (Subscript[f, 0])==(1-α) (Subscript[f, -1])+α (Subscript[f, 1]), where 0≤α≤1 is chosen so that (Subscript[f, 0]) is normal to the gradient ∇e, in particular α==((∇e·(Subscript[f, -1]))/(∇e·((Subscript[f, -1])-(Subscript[f, 1])))). 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 x'(t)==f(x(t)) with no explicit dependence on the independent variable t. The same ideas apply for a nonautonomous system x'(t)==f(t,x(t)) by adding a variable s(t)==t and considering the autonomous system {x'(t),s'(t)}=={f(s(t),x(t)),1}.

A simple example is provided by the vector field {1,sgn(1-y)} that has a discontinuity surface y==1. For this vector field (Subscript[f, 1])=={1,-1}, (Subscript[f, -1])=={1,1}, and e==1-y; therefore ∇e=={0,-1}, so α==(({0,-1}.{1,1})/({0,-1}.({1,1}-{1,-1})))==(1/2) and the sliding motion is governed by (Subscript[f, 0])==(1/2)(Subscript[f, -1])+(1/2)(Subscript[f, 1])=={1,0}.

Plot the vector field:

vf = VectorPlot[{1, Sign[1-y]},{x,0,2},{y,0,2}, VectorStyle->Gray, VectorPoints->2.Flatten[Outer[List, Range[0,11]/11,Range[0,11]/11],1]]

This attempts to solve the problem without any handling of the discontinuity surface:

fb[y_?NumberQ]:= Sign[1-y];
sol0 = First[NDSolve[{x'[t]==1,y'[t]==fb[y[t]],x[0]==0,y[0]==2},{x,y},{t,0,2}]];

ParametricPlot[Evaluate[{x[t],y[t]} /. sol0],{t,1,1.000036},AspectRatio->1]

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:

fb[y_?NumberQ]:= If[y > 1,-1,1];
sole = First[NDSolve[{x'[t]==1,y'[t]==fb[y[t]],x[0]==0,y[0]==2},{x,y},{t,0,2}, Method->"Extrapolation"]];

Show[vf, ParametricPlot[{x[t],y[t]} /. sole,{t,0,2}]]

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:

sol = First[NDSolve[{x'[t]==1,y'[t]==Sign[1-y[t]],x[0]==0,y[0]==2},{x,y},{t,0,2}]];

Show[vf, ParametricPlot[Evaluate[{x[t],y[t]} /. sol],{t,0,2}]]

When the discontinuous function is explicit in the equations, (RefLink[NDSolve, paclet:ref/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 smooth (ϵ)(z) that approaches sgn(z) as the parameter ϵ goes to zero:

smooth[ϵ_][z_]:= 2ArcTan[z/ϵ]/π
Plot[Evaluate[Table[smooth[(10^k)][z],{k,0,-4,-1/4}]],{z,-1,1}]

Manipulate the solution over a range of ϵ:

Manipulate[Module[{ϵ= 10.^k,ϵsol},
ϵsol = First[NDSolve[{x'[t]==1,y'[t]==smooth[ϵ][1-y[t]],x[0]==0,y[0]==2},{x,y},{t,0,2}]]; 
Show[vf,ParametricPlot[Evaluate[{x[t],y[t]} /. {sol,ϵsol}],{t,0,2}]]],{{k, -1, "(Subscript[Log, 10])[ϵ]"},-4,0}, SaveDefinitions->True]

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.

Exiting sliding mode.

In the figure above, the discontinuity surface is defined by e==0 and the vector fields for e>0 and e<0 are (Subscript[f, 1]) and (Subscript[f, -1]), respectively. Then on the surface ϕ use (Subscript[f, 0])==(1-α) (Subscript[f, -1])+α (Subscript[f, 1]), where 0≤α≤1 is chosen so that (Subscript[f, 0]) is normal to the gradient ∇e, in particular α==((∇e·(Subscript[f, -1]))/(∇e·((Subscript[f, -1])-(Subscript[f, 1])))). Beyond the point shown in both sides of the figure, ∇e.(Subscript[f, 1]) and ∇e.(Subscript[f, -1]) both have the same sign. On the left α>1 because ∇e.(Subscript[f, 1])>0 and ∇e.(Subscript[f, -1])>0, and on the right α<1 because ∇e.(Subscript[f, 1])<0 and ∇e.(Subscript[f, -1])<0, and in either case the trajectory exits sliding mode. These points can be found by numerically finding either the root α==1 or α==0 as appropriate.

Consider the vector field {1, ({{{, ({{y>1, -1}, {1-x/2, True}})}})}. Once x<2, 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:

sole = First[NDSolve[{x'[t]==1,y'[t]==If[y[t] - 1> 0, -1 , 1- x[t]/2],x[0]==0,y[0]==2},{x,y},{t,0,3}]];

Show[VectorPlot[{1, If[y >1, -1,1-x/2]},{x,0,3},{y,0,2}, VectorStyle->Gray, VectorPoints->Flatten[Outer[List, 3.Range[0,16]/16,2.Range[0,11]/11],1]], ParametricPlot[{x[t],y[t]} /. sole,{t,0,3}]]

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

Consider the vector field {y, ({{{, ({{-1, y>1}, {.1 y(t)-x(t), y<1}})}})}. For y<1, this is just an oscillator with negative damping. For y>1, all trajectories approach the discontinuity surface y==1. On y==1, 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:

Manipulate[Module[{x0, y0, sol},
{x0, y0} = p; sol =  First[NDSolve[{x'[t]==y[t],y'[t]==If[y[t] - 1>0, -1,.1 y[t] - x[t]],x[0]==x0,y[0]==y0},{x,y},{t,0,20}]];
ParametricPlot[{x[t], y[t]} /. sol,{t,0,20}, PlotRange->{{-1.5,1.5},{-1.5,1.5}}]],
{{p, {-1.5,1.5}},Locator}]

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 (((x+3/2))^2)+(y^2)+(((z-1))^2)==8 and the plane y==0 and the vector field given by {sgn(8-(((x+3/2))^2)-(y^2)-(((z-1))^2)),sgn(-y),1}.

This solves the system with two discontinuities:

sol = First[NDSolve[{x'[t]== Sign[8-y[t]^2-(x[t]+3/2)^2-(z[t]-1)^2],y'[t]==Sign[-y[t]],z'[t]==1,
x[0]==2,y[0]==-2, z[0]==0},{x,y,z},{t,0,5}]]
Plot[Evaluate[{x[t],y[t],z[t]}/.sol],{t,0,5}, ImageSize->Medium]

This shows the trajectory in 3D:

c1 = ContourPlot3D[8-y^2-(x+3/2)^2-(z-1)^2,{x,-2,2},{y,-2,2},{z,0,5}, Contours->0, Mesh->False];
c2 =  ContourPlot3D[-y,{x,-2,2},{y,-2,2},{z,0,5}, Contours->0,Mesh->False];

p1 = ParametricPlot3D[Evaluate[{x[t],y[t],z[t]} /. sol],{t,0,5}, PlotStyle->{Thickness[.015],Red}];
Show[c1,c2,p1, ImageSize->Medium]

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 (RefLink[NDSolve, paclet:ref/NDSolve]) does automatically for discontinuous functions.

Discontinuity signature variables are set up by including a (RefLink[WhenEvent, paclet:ref/WhenEvent]) expression where the event is on the discontinuity surface, and the action is rule setting the variable to "DiscontinuitySignature" and using the (RefLink[DiscreteVariables, paclet:ref/DiscreteVariables]) option to specify the range of the variable to be either {-1,0,1} or {-1,1}.

Setting up a discontinuity signature variable s[t].

Consider the first example in the previous section on sliding mode with the vector field {1,sgn(1-y)}.

Solve the system using an equivalent formulation with a discontinuity signature variable:

sol = First[NDSolve[{x'[t]==1,y'[t]==s[t],x[0]==0,y[0]==2,s[0]==-1,WhenEvent[1-y[t]==0,s[t]->"DiscontinuitySignature"]},{x,y,s},{t,0,2},DiscreteVariables->s[t]∈{-1,0,1}]]

Plot the solution and the values of the discontinuity signature variable s:

Row[{ParametricPlot[{x[t],y[t]} /. sol,{t,0,2}, PlotRange->{0,2}]Plot[s[t] /. sol, {t,0,2}, PlotStyle->Thickness[.02]]}]

The variable s switches from -1 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 d==0 and so sgn(d) is 0. In practice, the actual solution may not satisfy this due to numerical error. When in sliding mode, (RefLink[NDSolve, paclet:ref/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:

sol = First[NDSolve[{x'[t]==1,y'[t]==s[t],x[0]==0,y[0]==2,s[0]==-1,WhenEvent[1-y[t]==0,s[t]->"DiscontinuitySignature"]},{x,y,s},{t,0,2},DiscreteVariables->s[t]∈{-1,1}]]

The solution is incorrect:

ParametricPlot[{x[t],y[t]} /. sol,{t,0,2}]

An autonomous system defined by

can be represented by x'(t)==(Subscript[f, s(t)])(x(t)), where s(t) is set by the Wolfram Language expression (RefLink[WhenEvent, paclet:ref/WhenEvent])[e[x[t]]>0,s[t]->"DiscontinuitySignature"] in (RefLink[NDSolve, paclet:ref/NDSolve]).

The time derivative of the event function, ((de(x(t)))/(dt))==∇e.(Subscript[f, s])(x) is important in determining what happens when e(x((Subscript[t, a])))==0 at some time (Subscript[t, a]). There are six cases where s will be changed when e(x((Subscript[t, a])))==0 with (RefLink[WhenEvent, paclet:ref/WhenEvent])[e[x[t]]>0,s[t]->"DiscontinuitySignature"].

Crossing from below and above.

Case 1: (Cross from below) e<0 and s==-1 for t<(Subscript[t, a]) with ∇e.(Subscript[f, 1])(x)>0 at t==(Subscript[t, a]). The solution trajectory simply crosses to e>0 and s is set to 1.

Case 2: (Cross from above) e>0 and s==1 for t<(Subscript[t, a]) with ∇e.(Subscript[f, -1])(x)<0 at t==(Subscript[t, a]).The solution trajectory simply crosses to e<0 and s is set to -1.

Entering sliding mode from above and below.

Case 3: (Enter sliding mode from above) e>0 and s==1 for t<(Subscript[t, a]) with ∇e.(Subscript[f, -1])(x)>0 at t==(Subscript[t, a]). Since e is decreasing to zero it also holds that  ((de)/(dt))==∇e.(Subscript[f, 1])(x)≤ 0 at t==(Subscript[t, a]). The derivative ((de)/(dt)) 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 s is set to 0. The effective vector field (Subscript[f, 0])(x) is computed using the Filippov continuation (Subscript[f, 0])(x)==(Subscript[αf, 1])(x)+(1- α) (Subscript[f, -1])(x), where α==∇e.(Subscript[f, -1])(x)/(∇e.(Subscript[f, -1])(x)-∇e.(Subscript[f, 1])(x)). Note that 0≤α≤1 since ∇e.(Subscript[f, -1])(x)≥0 and ∇e.(Subscript[f, 1])(x)<0.

Case 4: (Enter sliding mode from below) e<0 and s==-1 for t<(Subscript[t, a]) with ∇e.(Subscript[f, 1])(x)<0 at t==(Subscript[t, a]). This is essentially the same as case 3 and the solution enters sliding mode with s set to 0.

Exiting sliding mode to above and below.

Case 5: (Exit sliding mode to above) e==0, s==0, ∇e.(Subscript[f, 1])(x)<0, and ∇e.(Subscript[f, -1])(x)>0 for t<(Subscript[t, a]) and at t==(Subscript[t, a]) ∇e.(Subscript[f, 1])(x)==0 with  (d/(dt))(∇e.(Subscript[f, 1])(x))> 0. The solution is no longer in sliding mode and s is set to 1. Note that the event that triggers this change is effectively a zero crossing of ∇e.(Subscript[f, 1])(x).

Case 6: (Exit sliding mode to below) e==0, s==0, ∇e.(Subscript[f, 1])(x)<0, and ∇e.(Subscript[f, -1])(x)>0 for t<(Subscript[t, a]) and at t==(Subscript[t, a]) ∇e.(Subscript[f, -1])(x)==0 with  (d/(dt))(∇e.(Subscript[f, 1])(x))< 0. The solution is no longer in sliding mode and s is set to -1. Note that the event that triggers this change is effectively a zero crossing of ∇e.(Subscript[f, -1])(x).

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

where (Subscript[f, 1]) and (Subscript[f, -1]) are constants that can be defined to see the different cases.

Define a function that gives the differential equations given definitions fa for (Subscript[f, 1]) and fb for (Subscript[f, -1]):

de[fa_, fb_] := {x'[t]==1,y'[t]==(1/2)(1+s[t])fa + (1/2)(1-s[t])fb};

Define a function solver that, given function definitions and an initial condition for y, solves the system using (RefLink[NDSolve, paclet:ref/NDSolve]) with s as a discontinuity signature variable:

solver[fa_,fb_][y0_]:= First[NDSolve[{de[fa, fb] ,x[0]==0,y[0]==y0,s[0]==0,WhenEvent[y[t]-1-Sin[π x[t]]/4==0,s[t]->"DiscontinuitySignature"]},{x,y,s},{t,0,2},DiscreteVariables->s[t]∈{-1,0,1}]];

Define a function that plots the solution superimposed on the vector field and a plot of the discontinuity signature variable s along with α if there is sliding mode:

plotter[fa_,fb_][y0_]:= Module[{sol, vf,psol, ps, sliding,g,α},
sol = solver[fa,fb][y0];
sliding = MemberQ[InterpolatingFunctionValuesOnGrid[s /. sol],0.];
psol = ParametricPlot[Evaluate[{{x[t],y[t]} /. sol,{t,1+Sin[π t]/4}}],{t,0,2},PlotStyle->{{Red, Thickness[.02]},{Black}},PlotRange->{{0,2},{0,2}}];
vf = VectorPlot[{1,If[y>1+Sin[π x]/4,fa,fb]},{x,0,2},{y,0,2},VectorStyle->Gray, VectorPoints->Flatten[Outer[List, 2.Range[0,16]/16,2.Range[0,11]/11],1],FrameLabel->{"x","y"},Epilog->Line[{{-1,1},{3,1}}]];
ps = Plot[s[t] /. sol,{t,0,2},PlotStyle->Thickness[.02],AxesLabel->{"t",If[sliding,"s,α","s"]},PlotRange->{-1,1}];
If[sliding,
g[x_,y_]=D[y-1-Sin[π x]/4,{{x,y}}];
α[s_ /; s==0,x_,y_] := ({1,fb}.g[x,y])/({0,fb-fa}.g[x,y]);
ps = Show[ps,Plot[α[s[t],x[t],y[t]] /. sol,{t,0,2},PlotStyle->Gray]]];
Row[{Show[vf,psol],"     ",ps}]]

Case 1: Cross from below with (Subscript[f, 1])==2 and (Subscript[f, -1])==1:

plotter[2,1][0]

Case 2: Cross from above with (Subscript[f, 1])==-1 and (Subscript[f, -1])==-2:

plotter[-1,-2][2]

Case 3: Enter sliding mode from above with (Subscript[f, 1])==-1 and (Subscript[f, -1])==1:

plotter[-1,1][2]

Case 4: Enter sliding mode from below with (Subscript[f, 1])==-1 and (Subscript[f, -1])==1:

plotter[-1,1][0]

Case 5: Exit sliding mode to above (Subscript[f, 1])==-1/4 and (Subscript[f, -1])==1:

plotter[-1/4,1][1]

Case 6: Exit sliding mode to below with (Subscript[f, 1])==-1 and (Subscript[f, -1])==1/3:

plotter[-1,1/3][2]

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:

sol=NDSolveValue[{y''[t]==-1+y'[t]^2,y[0]==1,y'[0]==0, WhenEvent[y[t]<0,{tend = t,"StopIntegration"}]},y,{t,0,Infinity}]

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

Show[Plot[sol[t],{t,0,tend},Frame->True,Axes->False,PlotStyle->Blue],
Graphics[{{Green,Circle[{0,1},0.025]},{Red,Circle[{tend,sol[tend]},0.025]}}]]

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, (RefLink[Automatic, paclet:ref/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 (Subscript[y, 2]) reaches its initial value and direction:

vsol = NDSolve[{(({{(Subscript[y, 1])'[t] == (Subscript[y, 2])[t] , (Subscript[y, 1])[0]==2}, {(Subscript[y, 2])'[t]==1000 (1 - (Subscript[y, 1])[t]^2) (Subscript[y, 2])[t] - (Subscript[y, 1])[t], (Subscript[y, 2])[0]==0}})),WhenEvent[(Subscript[y, 2])[t]<0,"StopIntegration"]},
{(Subscript[y, 1]), (Subscript[y, 2])},{t,3000}]

By selecting the endpoint of the (RefLink[NDSolve, paclet:ref/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 μ:

vper[μ_?NumberQ] := Module[{per = $Failed},
NDSolve[{(({{(Subscript[y, 1])'[t] == (Subscript[y, 2])[t] , (Subscript[y, 1])[0]==2}, {(Subscript[y, 2])'[t]==μ (1 - (Subscript[y, 1])[t]^2) (Subscript[y, 2])[t] - (Subscript[y, 1])[t], (Subscript[y, 2])[0]==0}})),
WhenEvent[(Subscript[y, 2])[t]<0,{per=t,"StopIntegration"}]},{},{t,Max[100, 3 μ]}];per]

This uses the function to compute the period at μ==1000:

vper[1000]

This shows a log-log plot of the period as a function of μ:

LogLogPlot[vper[μ],{μ,1,1000}]

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 (RefLinkPlain[EquationTrekker, paclet:EquationTrekker/tutorial/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:

system = GetNDSolveProblem["HenonHeiles"];
vars = system["DependentVariables"];eqns = {system["System"], system["InitialConditions"]}

The Poincaré section of interest in this case is the collection of points in the (Subscript[Y, 2])-(Subscript[Y, 4]) plane when the orbit passes through (Subscript[Y, 1])==0.

Since the actual result of the numerical integration is not required, it is possible to avoid storing all the data in (RefLink[InterpolatingFunction, paclet:ref/InterpolatingFunction]) by specifying the output variables list (in the second argument to (RefLink[NDSolve, paclet:ref/NDSolve])) as empty, or {}. This means that (RefLink[NDSolve, paclet:ref/NDSolve]) will produce no (RefLink[InterpolatingFunction, paclet:ref/InterpolatingFunction]) as output, avoiding storing a lot of unnecessary data. (RefLink[NDSolve, paclet:ref/NDSolve]) does give a message (RefLink[NDSolve, paclet:ref/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 (RefLink[WhenEvent, paclet:ref/WhenEvent]) that defines the Poincaré section. The event action is to use (RefLink[Sow, paclet:ref/Sow]) on the values of (Subscript[Y, 2]) and (Subscript[Y, 4]):

peqns = {eqns, WhenEvent[(Subscript[Y, 1])[T]==0,Sow[{(Subscript[Y, 2])[T], (Subscript[Y, 4])[T]}]]};

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 (RefLink[Sow, paclet:ref/Sow]) on the values of (Subscript[Y, 2]) and (Subscript[Y, 4]):

data=
Reap[
NDSolve[peqns,{}, {T, 10000},
Method->{"FixedStep","StepSize"->0.25,"Method"->{"ExplicitRungeKutta","DifferenceOrder"->4}},
 MaxSteps->∞];
];

This plots the Poincaré section. The collected data is found in the last part of the result of (RefLink[Reap, paclet:ref/Reap]), and the list of points is the first part of that:

psdata=data[[-1,1]];
ListPlot[psdata,Axes->False,Frame->True,AspectRatio->1]

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 (RefLink[Sow, paclet:ref/Sow]) on the values of (Subscript[Y, 2]) and (Subscript[Y, 4]):

sdata=
Reap[
NDSolve[peqns,{}, {T, 10000},
Method->{"FixedStep","StepSize"->0.25,"Method"->{"SymplecticPartitionedRungeKutta","DifferenceOrder"->4,"PositionVariables"->{(Subscript[Y, 1])[T], (Subscript[Y, 2])[T]}}}, MaxSteps->∞];
];

This plots the Poincaré section. The collected data is found in the last part of the result of (RefLink[Reap, paclet:ref/Reap]), and the list of points is the first part of that:

psdata=sdata[[-1,1]];
ListPlot[psdata,Axes->False,Frame->True,AspectRatio->1]

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:

system=GetNDSolveProblem["ArnoldBeltramiChildress"];
eqs = system["System"];
vars = system["DependentVariables"];
icvars = vars /. T->0;

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

Y1=eqs; Y1[[2,2]] = 0;Y1

Y2=eqs;Y2[[{1,3},2]] = 0; Y2

This defines the implicit midpoint method:

ImplicitMidpoint = {"ImplicitRungeKutta", "Coefficients"->"ImplicitRungeKuttaGaussCoefficients","DifferenceOrder"->2,"ImplicitSolver"->{"FixedPoint",AccuracyGoal->10, PrecisionGoal->10,"IterationSafetyFactor"->1}};

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

ABCSplitting = {"Splitting",
"DifferenceOrder"->2,
"Equations"->{Y2,Y1,Y2},
"Method"->{"LocallyExact",ImplicitMidpoint,"LocallyExact"}};

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

psect[ics_] := 
Module[{reapdata},
reapdata = 
Reap[
NDSolve[{eqs,Thread[icvars == ics], WhenEvent[(Subscript[Y, 2])[T]==0, Sow[{(Subscript[Y, 1])[T], (Subscript[Y, 3])[T]}]]},{},{T,1000},
 Method->{"FixedStep", "StepSize"->1/4,Method->ABCSplitting},MaxSteps->∞]
];
reapdata[[-1,1]]
];

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

data = 
Mod[Map[psect, {{4.267682454609692,0,0.9952906114885919},
{1.6790790859443243,0,2.1257099470901704},
{2.9189523719753327,0,4.939152797323216},
{3.1528896559036776,0,4.926744120488727},
{0.9829282640373566,0,1.7074633238173198},
{0.4090394012299985,0,4.170087631574883},
{6.090600411133905,0,2.3736566160602277},
{6.261716134007686,0,1.4987884558838156},
{1.005126683795467,0,1.3745418575363608},
{1.5880780704325377,0,1.3039536044289253},
{3.622408133554125,0,2.289597511313432},
{0.030948690635763183,0,4.306922133429981},
{5.906038850342371,0,5.000045498132029}}],
2π];
ListPlot[data,ImageSize->Medium]

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 (RefLink[NDSolve, paclet:ref/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 k of energy is left after a bounce:

Reflection[k_, ramp_][{x_, xp_, y_, yp_}] := Module[{gramp, gp, xpnew, ypnew},
gramp = -ramp'[x];
If[Not[NumberQ[gramp]],
Print["Could not compute derivative "];
Throw[$Failed]];
gramp = {-ramp'[x], 1};
If[ gramp.{xp,yp} == 0,
Print["No reflection"];
Throw[$Failed]];
gp= {1,-1} Reverse[gramp];
(k/(gramp.gramp)) (gp gp.{xp,yp} - gramp gramp.{xp, yp})]

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

BouncingBall[k_, ramp_,{x0_,y0_}] := 
Module[{data, end, bounces, prange, xmin,xmax,ymin,ymax,xtraj, ytraj, prev = {x0,y0}},
bounces = Part[
Reap[Sow[{x0,y0}];{xtraj, ytraj} = NDSolveValue[{x''[t]==0,x'[0]==0,x[0]==x0, y''[t]==-9.8, y'[0]==0,y[0]==y0,WhenEvent[y[t]<ramp[x[t]],{Sow[{x[t],y[t]}],If[Norm[prev - {x[t],y[t]}] < 10^-3 && Norm[{x'[t],y'[t]}] < 10^-3,end = t;"StopIntegration", prev = {x[t],y[t]}],{x'[t],y'[t]}->Reflection[k,ramp][{x[t],x'[t],y[t],y'[t]}]}]},{x,y},{t,0,10}]], 2, 1];
prange = {{xmin,xmax},{ymin,ymax}} = Outer[#2[bounces[[All,#1]]]&,{1,2}, {Min,Max}];
Show[Plot[ramp[x],{x,xmin, xmax},PlotRange->prange,Epilog->{PointSize[.025],Map[Point,bounces]}, AspectRatio->(ymax - ymin)/(xmax - xmin)], ParametricPlot[{xtraj[t],ytraj[t]},{t,0,end}, PlotStyle->Red]]]

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

ramp[x_] := If[x < 1,1 - x, 0];
BouncingBall[.7, ramp,{0,1.25}]

The ramp is now defined to be a quarter circle:

circle[x_] := If[x < 1, Sqrt[1 - x^2], 0];
BouncingBall[.7, circle,{.1,1.25}]

This adds a slight waviness to the ramp:

wavyramp[x_] := If[x < 1, 1 - x + .05 Cos[11 Pi x] ,0];
BouncingBall[.75, wavyramp,{0,1.25}]

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 (RefLink[WhenEvent, paclet:ref/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:

L = 50;
Timing[sgsol = First[NDSolve[{(∂_(t,t))u[t,x] ==(∂_(x,x))u[t,x] - Sin[u[t,x]],
 u[0,x] == (E^(-(((x-5))^2)))+(E^(-(((x+5))^2)/2)), Derivative[1,0][u][0,x] == 0, u[t,-L] == u[t,L],WhenEvent[Abs[u[t ,-L]]>0.01,"StopIntegration","LocationMethod"->"StepEnd"]}, u, {t,0,1000}, {x,-L,L}, Method->{"MethodOfLines",
"SpatialDiscretization"->{"TensorProductGrid", "DifferenceOrder"->"Pseudospectral"}}]]]

This extracts the ending time from the (RefLink[InterpolatingFunction, paclet:ref/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:

end = InterpolatingFunctionDomain[u /. sgsol][[1,-1]];
DensityPlot[u[t,x] /. sgsol,{x,-L,L}, {t,0, end}, Mesh->False, PlotPoints->100]

The "DiscretizedMonitorVariables" option affects the way the event is interpreted for PDEs; with the setting (RefLink[True, paclet:ref/True]), u[t,x] is replaced by a vector of discretized values. This is much more efficient because it avoids explicitly constructing the (RefLink[InterpolatingFunction, paclet:ref/InterpolatingFunction]) to evaluate the event:

L = 50;
Timing[sgsol = First[NDSolve[{(∂_(t,t))u[t,x] ==(∂_(x,x))u[t,x] - Sin[u[t,x]],
 u[0,x] == (E^(-(((x-5))^2)))+(E^(-(((x+5))^2)/2)), Derivative[1,0][u][0,x] == 0, u[t,-L] == u[t,L],WhenEvent[Abs[First[u[t ,x]]]>0.01,"StopIntegration","LocationMethod"->"StepEnd"]}, u, {t,0,1000}, {x,-L,L}, Method->{"MethodOfLines","DiscretizedMonitorVariables"->True,
"SpatialDiscretization"->{"TensorProductGrid", "DifferenceOrder"->"Pseudospectral"}}]]]

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:

μ=(1/82.45);
(SuperStar[μ])=1-μ;
(Subscript[r, 1])=(Sqrt[((((Subscript[y, 1])[t]+μ))^2)+(((Subscript[y, 2])[t])^2)]);
(Subscript[r, 2])=(Sqrt[((((Subscript[y, 1])[t]-(SuperStar[μ])))^2)+(((Subscript[y, 2])[t])^2)]);
eqns={{(Subscript[y, 1])'[t]==(Subscript[y, 3])[t],(Subscript[y, 1])[0]==1.2},{(Subscript[y, 2])'[t]==(Subscript[y, 4])[t],(Subscript[y, 2])[0]==0},{(Subscript[y, 3])'[t]==2 (Subscript[y, 4])[t]+(Subscript[y, 1])[t]-(((SuperStar[μ]) ((Subscript[y, 1])[t]+μ))/(Subscript[r, 1]^3))-((μ ((Subscript[y, 1])[t]-(SuperStar[μ])))/(Subscript[r, 2]^3)),(Subscript[y, 3])[0]==0},{(Subscript[y, 4])'[t]==-2 (Subscript[y, 3])[t]+(Subscript[y, 2])[t]-(((SuperStar[μ]) (Subscript[y, 2])[t])/(Subscript[r, 1]^3))-((μ (Subscript[y, 2])[t])/(Subscript[r, 2]^3)),(Subscript[y, 4])[0]==-1.0493575098303198}};

The event function is the derivative of the distance from the initial conditions. A local maximum or minimum occurs when the value crosses zero:

ddist = 2((Subscript[y, 3])[t]((Subscript[y, 1])[t]-1.2)+(Subscript[y, 4])[t](Subscript[y, 2])[t]);

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 (RefLink[Evaluate, paclet:ref/Evaluate]) is necessary in the first argument of (RefLink[WhenEvent, paclet:ref/WhenEvent]) so that the variable substitutions get done correctly:

sol = First[NDSolve[{eqns,WhenEvent[Evaluate[ddist<0], tfar = t],WhenEvent[Evaluate[ddist>0],{tfinal=t,"StopIntegration"}]},{(Subscript[y, 1]),(Subscript[y, 2]),(Subscript[y, 3]),(Subscript[y, 4])},{t,∞}]]

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:

ParametricPlot[{(Subscript[y, 1])[t], (Subscript[y, 2])[t]} /. sol, {t,0,tfinal}, Epilog->{PointSize[.025],Point[{(Subscript[y, 1])[tfar], (Subscript[y, 2])[tfar]} /. sol]}]

Infection Model with Delay

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

system = {y1'[t]==-y1[t] y2[t-1]+y2[t-10],y1[t/;t≤0]==5,y2'[t]==y1[t] y2[t-1]-y2[t],y2[t/;t≤0]==1/10,y3'[t]==y2[t]-y2[t-10],y3[t/;t≤0]==1};
vars = {y1[t],y2[t],y3[t]};

Define a function that makes a (RefLink[WhenEvent, paclet:ref/WhenEvent]) expression to capture the local maximum of the i(^th) component. A separate tag for (RefLink[Sow, paclet:ref/Sow]) is used to distinguish the components:

lmax[i_] := With[{vi = vars[[i]], vpi = D[vars[[i]],t]},WhenEvent[vpi<0,Sow[{t,vi},i]]];

Solve the system with a (RefLink[WhenEvent, paclet:ref/WhenEvent]) for the maximum of each component:

data = 
Reap[
sol = First[NDSolve[{system,Array[lmax,3]},vars,{t,0,40}]],
{1,2,3}
]

Display the local maxima together with the solution components:

colors = {{Red},{Blue},{Green}};
plots = Plot[Evaluate[vars /. sol],{t,0,40},PlotStyle->colors];
max = ListPlot[Part[data,-1,All,1],PlotStyle->colors];
Show[plots,max]

Discontinuity Signature and C Code

For symbolically specified differential equations, (RefLink[NDSolve, paclet:ref/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 (RefLink[LibraryFunction, paclet:ref/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:

csource = "
/* Include required header */
#include \"WolframLibrary.h\"
#include \"math.h\"
#include \"stdlib.h\"

/* Return the version of Library Link */
DLLEXPORT mint WolframLibrary_getVersion() {
	return WolframLibraryVersion;
}

DLLEXPORT int slidingrhs(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res)
{
	int err;
	mint p, s;
	mreal t, x, y, fx, fy;
	MTensor Txy, Trhs;
	if (Argc != 3) return LIBRARY_FUNCTION_ERROR;
	t = MArgument_getReal(Args[0]);
	Txy = MArgument_getMTensor(Args[1]);
	s = MArgument_getInteger(Args[2]);
	if (libData->MTensor_getType(Txy) != MType_Real) 
		return LIBRARY_TYPE_ERROR;
	if (libData->MTensor_getRank(Txy) != 1) return LIBRARY_RANK_ERROR;
	p = 1;
	err = libData->MTensor_getReal(Txy, &p, &x);
	if (err) return err;
	p = 2;
	err = libData->MTensor_getReal(Txy, &p, &y);
	if (err) return err;
	
	fx = x + y;
	if (s == 1) {
		fy = -2.*t;
	} else if (s == -1) {
		fy = 3.*t;
	} else {
		return LIBRARY_NUMERICAL_ERROR;
	}
	    
	p = 2;
	err = libData->MTensor_new(MType_Real, 1, &p, &Trhs);
	if (err) return err;

	p = 1;
	err = libData->MTensor_setReal(Trhs, &p, fx);
	if (err) return err;

	p = 2;
	err = libData->MTensor_setReal(Trhs, &p, fy);
	if (err) return err;

	MArgument_setMTensor(Res, Trhs);
    return 0;
}
";

Use (PackageLink[CreateLibrary, paclet:CCompilerDriver/ref/CreateLibrary, CCompilerDriver Package Symbol]) from the CCompilerDriver package to make a library:

Needs["CCompilerDriver`"]

lib = CreateLibrary[csource, "sliding"]

Load the library into a (RefLink[LibraryFunction, paclet:ref/LibraryFunction]):

lfun = LibraryFunctionLoad[lib, "slidingrhs", {Real, {Real, 1}, Integer}, {Real, 1}]

Solve the system for X=={x,y} with the switching on the circle (x^2)+(y^2)==1:

sol = First[NDSolve[{X'[t] == lfun[t, X[t], s[t]], X[0]=={0,0}, s[0]== -1, WhenEvent[X[t].X[t]- 1==0, s[t]->"DiscontinuitySignature"]}, {X, s}, {t,0,2}, DiscreteVariables->{Element[s[t], {-1,0,1}]}]]

Show a plot of the solution with blue for s<0, green for s>0, and red for sliding mode:

plot = Show[RegionPlot[x^2 + y^2 ≤ 1, {x,0,1},{y,0,1}],ParametricPlot[Evaluate[If[s[t]<0, X[t]] /. sol], {t,0,2}, PlotStyle->{{Blue, Thickness[0.015]}}],ParametricPlot[Evaluate[If[s[t]==0, X[t]] /. sol], {t,0,2}, PlotStyle->{{Red, Thickness[0.015]}}],ParametricPlot[Evaluate[If[s[t]>0, X[t]] /. sol], {t,0,2}, PlotStyle->{{Green, Thickness[0.015]}}], AspectRatio->Automatic]

With automatic discontinuity handling, an equivalent solution can be found with symbolic input:

ssol = First[NDSolve[{x'[t]==x[t] + y[t] , y'[t]==If[(x[t]^2 + y[t]^2)> 1, -2 t, 3 t], x[0]==0, y[0]==0}, {x,y}, {t,0,2}]];
Show[plot, ParametricPlot[{x[t], y[t]} /. ssol, {t,0,2}, PlotStyle->White]]