Solution of Differential Equations with Discontinuities

Rob Knapp

Wolfram Research, Inc.

Introduction

The methods for NDSolve typically assume some level of smoothness of the right hand side of the differential quations.  When this assumption is violated, the methods can be inefficient, or even give very inaccurate results.

A method for handling discontinuities by treating them as events is currently under development.  The method is designed to use symbolic analysis to automatically detect both impulse (DiracDelta) and jump discontinuities in the differential equation.  In this presentation, I will show you several examples of the new method at work.

Stick-slip Example

Drag a block (position x(t)) across a surface by moving the anchor point of a spring (position a(t)).
Assume that static friction is exceeds the sliding friction at low velocities.

Starting from rest, initially the block will remain stuck until the force from the spring is sufficient to break the static friction, at which point it begins to slip.

If the movement of the spring anchor is slow enough, the action of the spring will be such that the velocity is zero again at some point.  Then, the static friction takes over again until the spring force is sufficient to break it.

μ = .25; (* Coefficient of sliding friction *)
 α = 0.001; 
(* Spring force S -- Nonlinearity to prevent block from passing the anchor *)
(Subscript[F, spring])[a_, x_]:= (a - x - 1)(1 + 1/(a - x)); 
(* Friction *)
(Subscript[F, static])[] :=2 μ; 
(Subscript[F, sliding])[xdot_] := -μ xdot;

The following function gives the net force on the block.

If the block is moving (above a threshold speed), then the force is given by the spring force minus the sliding friction.

If the block is static, the static friction will be equal and opposite to the spring force up to the maximum of the static friction force.

Force[a_,x_,xdot_] := EvaluatedIf[(xdot^2 - α^2) > 0, 
(* Sufficient velocity for sliding friction *)
(Subscript[F, spring])[a,x] + (Subscript[F, sliding])[xdot],
(* else *)
EvaluatedIf[((Subscript[F, spring])[a,x]^2 - (Subscript[F, static])[]^2) > 0,
(* Sufficient force to overcome static friction *)
(Subscript[F, spring])[a,x] + (Subscript[F, sliding])[xdot],
(* else static: damp any residual velocity quickly *)
-10xdot]];

Here is what heppens with the default (LSODA) method

sol = First[NDSolve[{
x''[t] == Force[a[t],x[t],x'[t]],x[0]== x'[0]==0,
a'[t] == .1, a[0]== 1},
{x,a},{t,0,30}]]

Length[times = StepTimes[sol]]

Show[Plot[Evaluate[{Tooltip[x[t],"x"],Tooltip[x'[t],"x'"],Tooltip[ x''[t],"x''"],Tooltip[(Subscript[F, spring])[a[t], x[t]]^2 - 4 μ^2,"Static Friction Threshold"],Tooltip[(x'[t])^2 - α^2,"Veocity Threshold"] } /. sol], {t,0,30}, PlotRange->{-1,1}],ListPlot[Transpose[{times, x''[times] /. sol}]]]

This allows zooming in:

DynamicModule[{plotrange={{0,30},{-1,1}},pt={0,0},beginPt,rubberband={}},LocatorPane[Dynamic[pt,{(beginPt=#)&,(rubberband=Rectangle[beginPt,#])&,(rubberband={};plotrange=Transpose[Sort[{beginPt,#}]])&}],Dynamic[
{tmin, tmax} = plotrange[[1]];
tl = Select[times, Function[tmin ≤ # ≤ tmax]];Show[Plot[Evaluate[{Tooltip[x[t],"x"],Tooltip[x'[t],"x'"],Tooltip[ x''[t],"x''"],Tooltip[(Subscript[F, spring])[a[t], x[t]]^2 - 4 μ^2,"Static Friction Threshold"],Tooltip[(x'[t])^2 - α^2,"Veocity Threshold"] } /. sol], {t,tmin, tmax}, PlotRange->plotrange],ListPlot[Transpose[{times, x''[times] /. sol}]],Epilog->{FaceForm[],EdgeForm[Thin],rubberband}]],Appearance->None]]

And the new method:

dsol = First[NDSolve[{
x''[t] == Force[a[t],x[t],x'[t]],x[0]== x'[0]==0,
a'[t] == .1, a[0]== 1},
{x,a},{t,0,30},Method->"DiscontinuityHandler"]]

Length[dtimes = StepTimes[dsol]]

Show[Plot[Evaluate[{Tooltip[x[t],"x"],Tooltip[x'[t],"x'"],Tooltip[ x''[t],"x''"],Tooltip[(Subscript[F, spring])[a[t], x[t]]^2 - 4 μ^2,"Static Friction Threshold"],Tooltip[(x'[t])^2 - α^2,"Veocity Threshold"] } /. dsol], {t,0,30}, PlotRange->{-1,1}],ListPlot[Transpose[{times, x''[times] /. dsol}]]]

This allows zooming in:

DynamicModule[{plotrange={{0,30},{-1,1}},pt={0,0},beginPt,rubberband={}},LocatorPane[Dynamic[pt,{(beginPt=#)&,(rubberband=Rectangle[beginPt,#])&,(rubberband={};plotrange=Transpose[Sort[{beginPt,#}]])&}],Dynamic[
{tmin, tmax} = plotrange[[1]];
tl = Select[times, Function[tmin ≤ # ≤ tmax]];Show[Plot[Evaluate[{Tooltip[x[t],"x"],Tooltip[x'[t],"x'"],Tooltip[ x''[t],"x''"],Tooltip[(Subscript[F, spring])[a[t], x[t]]^2 - 4 μ^2,"Static Friction Threshold"],Tooltip[(x'[t])^2 - α^2,"Veocity Threshold"] } /. dsol], {t,tmin, tmax}, PlotRange->plotrange],ListPlot[Transpose[{dtimes, x''[dtimes] /.dsol}]],Epilog->{FaceForm[],EdgeForm[Thin],rubberband}]],Appearance->None]]

This makes an animation of the stick-slip phenomenon

b = Table[Blocks[x[t] /. dsol, a[t] /. dsol, PlotRange->{{0,4},{0,1}}],{t,0,30,.1}];
ListAnimate[b]

Another interesting case is when the anchor is moved periodically

psol = First[NDSolve[{
x''[t] == Force[a[t],x[t],x'[t]],x[0]== x'[0]==0,
a'[t] ==  .25 Cos[.25 t], a[0]==1},
{x,a},{t,0,100}, Method->"DiscontinuityHandler"]]

StepPlot[psol,{x,x',x''}]

b = Table[Blocks[x[t] /. psol, a[t] /. psol,PlotRange->{{-2,2},{0,1}}],{t,0,100,.25}];
ListAnimate[b]

Impulse Discontinuities

I will use "impulse discontinuity" to refer to a situation where the solution itself has a jump at a certain time.  Typically, the impulse can be expressed by a dirac delta function.

sol = First[NDSolve[{x''[t] + Sin[x[t]] == 0.1 DiracDelta[Cos[t/4]],x[0]==x'[0]==0},x,{t,0,20}]]

The default method just uses sampling, so the delta function is entirely missed.

StepPlot[sol,{x,x',x''}]

With the new method, the presence of delta functions is detected symbolically, with each instance being set up as an event to locate:

dsol = First[NDSolve[{x''[t] + Sin[x[t]] == 0.1DiracDelta[Cos[t/4]],x[0]==x'[0]==0},x,{t,0,1000},Method->"DiscontinuityHandler"]]

StepPlot[dsol,{x,x',x''},{0,20}]

StepPhasePlot[dsol,{x,x'}]

Bouncing Ball

This appears in the advanced documentation for the NDSolve method EventLocator as an example of how you can use events to stop and restart the integration.  With the DiscontinuityHandler event, there is an easier, more intuitive way to do this.

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];
{xpnew, ypnew} = (k/(gramp.gramp)) (gp gp.{xp,yp} - gramp gramp.{xp, yp});
{x, xpnew, y, ypnew}]

Bounce[{k_, ramp_},t_?NumberQ,{x_,xp_,y_,yp_}]:= If[1-tprev/t<0.01,Throw[Null,"StopIntegration"],tprev=t;(Reflection[k,ramp][{x,xp,y,yp}]-{x,xp,y,yp})]

MakeBounces[ramp_, k_,{x0_,y0_}] := Block[{tprev = 0, t, x, y},
bsol = First[NDSolve[{x''[t]==0,y''[t]==-9.8,x[0]==x0,x'[0]==0,y[0]==y0,y'[0]==0},{x,y},{t,0,10},Method->{"DiscontinuityHandler","Discontinuities":>y[t]-ramp[x[t]],"DiscontinuityAction"->Bounce[{k,ramp},t,{x[t],x'[t],y[t],y'[t]}]}, MaxStepSize->0.01]];
end = Part[x /. bsol, 3,1,-1];pp = Show[ParametricPlot[Evaluate[{x[t],y[t]} /. bsol],{t,0,end},PlotRange->All],Plot[Evaluate[ramp[x]],{x,0,1},PlotStyle->Green, Axes->False]]];

This sets up a simple interface for bouncing the ball

DynamicModule[{ramp, k = 0.75, at = 0, wind,pt = {0,1.25}},
end = 1;
r[1][x_] := If[x < 1,1 - x, 0];
r[2][x_] := If[x < 1, Sqrt[1 - x^2], 0];
r[3][x_] := If[x < 1, 1 - x + .05 Cos[11 Pi x] ,0];
Column[{TabView[Table[{r[i],Plot[r[i][x],{x,0,1},AspectRatio->1,ImageSize->Tiny,Axes->False]->r[i][x]},{i,3}],Dynamic[ramp,Function[Set[ramp,#];MakeBounces[ramp, k, pt]]]],Slider[Dynamic[k,Function[Set[k,#];MakeBounces[ramp, k, pt]]]], Animator[Dynamic[at]],Dynamic[Show[pp,Graphics[{Locator[Dynamic[pt,Function[Set[pt,#];MakeBounces[ramp, k, pt]]]], {Red,PointSize[0.025],Point[{x[end at],y[end at]} /. bsol]}}],PlotRange->{{-.1,2},{1.3,0}}, Frame->True]]}]]

Jump Discontinuities

I will use "Jump Discontinuity" for cases where the rhs of the DE's have a jump, so the solutions are continuous, but the derivatives are not necessarily so.

sol = First[NDSolve[{x''[t] + 0.1 x'[t] + Round[3 Sin[x[t]]]/3==Cos[t],x[0]== 0,x'[0]==0},x,{t,0,100}]]

StepPlot[sol,x'',{0,20}]

StepPhasePlot[sol,{x,x'}]

With the new method, it is possible to integrate out much farther in time.

dsol = First[NDSolve[{x''[t] + 0.1 x'[t] + Round[3 Sin[x[t]]]/3==Cos[t],x[0]== 0,x'[0]==0},x,{t,0,2000},Method->"DiscontinuityHandler",MaxSteps->100000]]

StepPlot[dsol,x'',{0,20}]

The solution eventually settles down towards some sort of orbit

StepPhasePlot[dsol, {x,x'},{1900, 2000}]

StepPlot[dsol,x'', {1950, 2000}]

Sliding Motion

When a discontinuity surface is encountered where the vector field on both sides of the discontinuity points towards the discontinuity, the state can enter  "sliding motion" where the solution effectively slides along the discontinuity surface.  This was fir defined by Filippov and is sometimes referred to as Filippov sliding motion.  Filippov's big contribution was to observe that the when the condition that the vector field points in on both sides is met, the solution can be continued by taking a convex combination of the vector fields on either side that lies in the surface.

Manipulate[
{x1,y1} = p1;
{x2,y2} = p2;
p3 = If[y1 > 0,p1,If[y2 < 0,p2, λ =y2/(y2 -y1); λ p1 + (1 - λ)p2]];
Graphics[{{Green, Tooltip[Arrow[{{0,0}, p1}],"Above"]},{Blue, Tooltip[Arrow[{{0,0},p2}],"Below"]}, {Red, Thickness[0.012],Arrow[{{0,0}, p3}]},{Black, Line[{{-1.1,0},{1.1,0}}]}},PlotRange->{{-1.1,1.1},{-1.1,1.1}}],
{{p1,{1,-1}},Locator}, {{p2,{1,1}}, Locator}]

Sliding Motion Example

Consider an example where above the line y = 1, the vector field points down and to the right and below the line y = 1 the vector field is from negatively damped harmonic motion.  At several points along the discontinuity line, the conditions for sliding motion are satisfied.

Show[vf = DrawVectorField2D[{y,If[y > 1, -1,.1 y - x]},{{x,-2,1.1},{y, -1.25,1.5}},ScalingFactor->0.05],Frame->True]

The default method gets stuck:

sol = First[NDSolve[{x'[t]==y[t],y'[t]==If[y[t] - 1>0, -1, .1 y[t] - x[t]],x[0]==-2,y[0]==1.5},{x,y},{t,0,20}]]

The solution is just chattering around the discontinuity.

Show[StepPhasePlot[sol, {x,y},{.5,Max}],AspectRatio->1]

The new method has no problem with the solution:

dsol = First[NDSolve[{x'[t]==y[t],y'[t]==If[y[t] - 1>0, -1, .1 y[t] - x[t]],x[0]==-2,y[0]==1.5},{x,y},{t,0,20}, Method->"DiscontinuityHandler"]]

StepPhasePlot[dsol,{x,y}]

In some ways, the sliding surface acts like a very strong attractor

Manipulate[dsol = First[NDSolve[{x'[t]==y[t],y'[t]==If[y[t] - 1>0, -1, -γ y[t] - x[t]],x[0]==p[[1]],y[0]==p[[2]]},{x,y},{t,0,20}, Method->"DiscontinuityHandler"]];
StepPhasePlot[dsol, {x,y}, All, PlotRange->{{-3,2},{-2,2}}],
{{p, {-2,1.5}},Locator},{{γ,-.1},-1,1}]

Pseudo Sliding Motion

You are quite likely wondering if this continued solution is really justified.  One way to see this is to consider it as the limit of a sequence of solutions with successively sharper discontinuities.  Pseudo sliding motion is what you get when you have a smooth, but steep change with the vector fields pointing in across the surface of change.

Consider the smoothed unit step function,

(Subscript[s, ϵ_])[x_] := (ArcTan[x/ϵ]/π + 1/2)

You can see that it approaches a unit step function as ϵ-> 0

Manipulate[Plot[{UnitStep[x],(Subscript[s, ϵ])[x]},{x,-1,1}, Exclusions->None],{{ϵ,1},1.*^-4,1}]

Now compare sliding mode versus the solution with the unit step function replaced by (Subscript[s, ϵ]).

imethod = {"ExplicitRungeKutta", "DifferenceOrder"->8, "StiffnessTest"->False};

sol = First[NDSolve[{x'[t]==y[t],y'[t]==UnitStep[1 - y[t]](.1 y[t] - x[t] + 1) - 1,x[0]==-2,y[0]==1.5},{x,y},{t,0,10},Method->{"DiscontinuityHandler", Method->imethod}]];

Moving the slider allows you to see how the solution approaches sliding mode.

{{xmin,xmax},{ymin,ymax}} = {{-2,1.5},{-1.5,1.5}};
Manipulate[
psol = First[NDSolve[{x'[t]==y[t],y'[t]==(Subscript[s, ϵ])[1 - y[t]](.1y[t]-x[t] + 1) - 1,x[0]==-2,y[0]==1.5},{x,y},{t,0,10},Method->imethod,InterpolationOrder->All, MaxStepSize->100/ϵ]];
vf = DrawVectorField2D[{y,(Subscript[s, ϵ])[1 - y](.1y-x + 1) - 1},{{x,xmin, xmax},{y, ymin, ymax}}];
pl = ParametricPlot[Evaluate[{x[t],y[t]} /.{sol, psol}],{t,0,10}, PlotRange->{{xmin,xmax},{ymin,ymax}}];
Show[pl,vf],{{ϵ,1},1.*^-4,1}]

A close up helps show how the vector field is changing ...

{{xmin,xmax},{ymin,ymax}} = {{-2,1.5},{.9,1.1}};
Manipulate[
psol = First[NDSolve[{x'[t]==y[t],y'[t]==(Subscript[s, ϵ])[1 - y[t]](.1y[t]-x[t] + 1) - 1,x[0]==-2,y[0]==1.5},{x,y},{t,0,10},Method->imethod,InterpolationOrder->All, MaxStepSize->100/ϵ]];
vf = DrawVectorField2D[{y,(Subscript[s, ϵ])[1 - y](.1y-x + 1) - 1},{{x,xmin, xmax},{y, ymin, ymax}}];
pl = ParametricPlot[Evaluate[{x[t],y[t]} /.{sol, psol}],{t,0,10}, PlotRange->{{xmin,xmax},{ymin,ymax}}, AspectRatio->1];
Show[pl,vf],{{ϵ,1},1.*^-4,1}]

Summary

A new method is under development for handling differential equations with discontinuities in the ride hand side.

It processes differential equations symbolically to discover

Impulse discontinuities

Jump discontinuities

When the jump leads to sliding motion, the appropriate convex combination is computed automatically.

Though the default method can sometimes get a solution, the new method will typically get a more accurate solution in far fewer steps.

Initialization

EvaluatedIf[args___] := If[args]; (* Used so that derivatives are computed symbolically *)

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]

StepTimes[{_->if_InterpolatingFunction,___Rule}] := First[InterpolatingFunctionCoordinates[if]];

TimeRange[{tmin_, tmax_}, times_] := Module[{trmin = tmin, trmax = tmax},
If[trmin == Min, trmin = times[[1]]];
If[trmax == Max, trmax = times[[-1]]];
{trmin, trmax}];
TimeRange[All, times_] := times[[{1,-1}]];

StepPlot[sol:{(_->_InterpolatingFunction)..}, v_List, tr_, opts___] := Module[{t, times = StepTimes[sol],trmin, trmax},
{trmin, trmax} = TimeRange[tr, times];
times = Select[times, Function[trmin ≤ # ≤ trmax]];
Show[Plot[Evaluate[Map[Tooltip[#[t],ToString[#]]&,v] /. sol],{t,trmin, trmax}, opts], ListPlot[Map[Transpose[{times, #[times] /. sol}]&, v]]]];
StepPlot[sol:{(_->_InterpolatingFunction)..}, v_List]  := StepPlot[sol, v, All];
StepPlot[sol:{(_->_InterpolatingFunction)..}, v_, rest___]:= StepPlot[sol, {v}, rest];

StepPhasePlot[sol:{(_->_InterpolatingFunction)..},{x_,y_}, tr_, opts___] :=Module[{t, times = StepTimes[sol],trmin, trmax},
{trmin, trmax} = TimeRange[tr, times];
times = Select[times, Function[trmin ≤ # ≤ trmax]];
Show[ParametricPlot[Evaluate[{x[t],y[t]} /. sol],{t,trmin,trmax},opts], Graphics[{Red, Point[Transpose[{x[times],y[times]} /. sol]]}]]];
StepPhasePlot[sol:{(_->_InterpolatingFunction)..},{x_,y_}] := StepPhasePlot[sol,{x,y},All];

Blocks[x_, a_, opts___] := Module[{s = .25},Graphics[{{Blue, Thickness[.01/Abs[x - a]], Line[{{x,s},{a,s}}]},{Red, Polygon[{{x - s,0},{x+s,0},{x+s,2 s},{x-s,2 s}}]},{Cyan,Disk[{a,s},s/2]}},opts, Axes->True]];

Options[DrawVectorField2D] = {PlotPoints->24, ScalingFactor->Automatic};

DrawVectorField2D[fofxy_, {{x_, xmin_, xmax_},{y_, ymin_, ymax_}},OptionsPattern[]] := Module[{f, pp, ppx, ppy, sc, X, Y, pxy},
f = Experimental`CreateNumericalFunction[{x, y},fofxy,{2}];
pp = OptionValue[PlotPoints];
If[Not[ListQ[pp]],pp = {pp, pp}];
If[Length[pp] ≠ 2, Return[$Failed]];
If[Not[VectorQ[pp, Function[Positive[#] && IntegerQ[#]]] && (Length[pp]==2)], 
Message[MessageName[DrawVectorField2D,"plotpoints"], pp];
Return[$Failed]
];
{ppx, ppy} = pp;
sc = OptionValue[ScalingFactor];
dx = (xmax - xmin)/ppx;
X = xmin + dx(Range[0, ppx-1] + .5);
dy = (ymax-ymin)/ppy;
Y = ymin + dy(Range[0, ppy-1] + .5);
fxy = Flatten[Outer[f[{##}]&, X, Y], 1];
If[sc === Automatic, 
sx = dx/Max[Abs[fxy[[All,1]]]];
sy = dy/ Max[Abs[fxy[[All,2]]]];
sc = .5 Min[sx, sy]
];
xy = Flatten[Outer[List, X, Y],1];
fxy = xy + sc*fxy;
Graphics[{Arrowheads[Tiny],MapThread[Arrow[{##}]&, {xy, fxy}]}]
]