Event Handling in Differential Equations

The following animation shows a ball bouncing down a nonlinear slope (a hybrid system). When the ball hits the slope, NDSolve automatically stops the computation, resets the conditions for bouncing off, and restarts. Detected events have occurred at the gray spots in the plot.

BouncingBallAnimation[
        If[#  <  1,  1  -  #  +    .05 Cos[11 Pi #]  , 0] &, {0, 
          1.25},  {0, 0},  DragCoefficient -> 1,  
        RestitutionCoefficient ->  .9, 
        Fills -> {RGBColor[1,  .3, 0]}];

Initialization

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]

Needs["Graphics`FilledPlot`"];

\(\(OneBounce[surf_, \ k_, \ ν_]\)[{t0_, x0_, \ 
          xp0_}]\  := \ 
Module[{sol, \ t1, \ x1, \ x1p, \ 
          g}, 
g[t_?NumberQ]\  = 
          Module[{gv\  = \ 0. \ x0}, \ gv[\([\(-1\)]\)]\  = \ \(-9.8\); \ 
            gv]; 
sol\  = \ 
          First[NDSolve[{\(x''\)[t]\  == 
                  g[t]\  - \ ν\ \ Norm[\(x'\)[t]]\ \(x'\)[t]\ , \ 
                x[t0]\  == \ x0, \ \(x'\)[t0]\  == \ xp0}, 
              x, {t, t0, \ ∞}, \ 
              Method -> {EventLocator, \ "Event" :> \ 
                    Last[x[t]]\  - \ 
                      Apply[surf, \ 
                        Most[x[t]]]}, 
MaxStepSize -> 
                0.01]]; 
t1\  = \ \(InterpolatingFunctionDomain[
              x\  /. \ sol]\)[\([1, \(-1\)]\)]; 
{x1, \ 
            xp1}\  = \ {x[t1], \ \(x'\)[t1]}\  /. \ 
            sol; 
xp1 = \ \(Reflection[surf, \ k]\)[{x1, \ 
              xp1}]; 
Sow[{x[#]\  /. sol, 
            t0\  ≤ \ #\  ≤ \ \ t1}, \ "X"]; 
Sow[x1, \ "Bounces"]; 
{t1, \ x1, 
          xp1}]\)

\(\(Reflection[surf_, \ k_]\)[{x_, \ xp_}]\  := \ 
      Module[{zeroes, \ gsurf, \ normal, \ xpnew, \ 
          n\  = \ Length[x]\  - \ 1}, 
zeroes\  = \ 
          Table[0, \ {n}]; 
gsurf\  = \ 
          Table[
der\  = \ Apply[Derivative, zeroes]; 
            der[\([i]\)]\  = \ 1; 
Apply[der[surf], \ 
              Most[x]], 
{i, n}]; 
If[
          Not[VectorQ[gsurf, \ 
              NumberQ]], 
Print["Could not compute derivative ", \ 
            gsurf]; 
Throw[$Failed]]; 
normal\  = \ 
          Append[\(-gsurf\), \ 1]; 
normal\  = \ 
          normal/Norm[normal]; 
If[\ 
          normal . xp\  == \ 
            0, 
Print["No reflection"]; 
Throw[$Failed]]; 
xpnew\  = \ 
          k\ \((xp\  - \ 
                2\ normal\ normal . xp)\)
]\)

Options[BouncingBall]  =  {DragCoefficient ->  .1, 
          RestitutionCoefficient ->  .9};

\(BouncingBallAnimation[surf_, \ x0_, \ xp0_, \ 
        opts___]\  := \ 
Module[{flops, \ ν, \ k, \ 
          data, bounces, \ solution, \ nw, 
          result}, 
flops\  = \ 
          Flatten[{opts, \ 
              Options[BouncingBall]}]; 
ν\  = \ 
          N[DragCoefficient\  /. \ 
              flops]; 
If[\(! TrueQ[ν\  ≥ \ 0]\), \ 
Print["Drag coefficient must be nonnegative"]; Return[$Failed]]; 
k\  = \ 
          N[RestitutionCoefficient\  /. \ 
              flops]; 
If[\(! TrueQ[
              0\  < \ k\  ≤ \ 
                1]\), \ 
Print["Restitution coefficient must a positive number <= 1"]; Return[$Failed]]; 
If[
          Apply[surf, \ Most[x0]]\  ≥ \ 
            Last[x0], 
Print["Start above the ramp"]; 
Return[$Failed]]; 
data\  = 
          Reap[
Catch[
              Sow[x0, \ "Bounces"]; 
nw\  = \ 
                NestWhile[OneBounce[surf, \ k, \ ν], \ {0, x0, \ xp0}, \ 
                  Function[1\  - \ #1[\([1]\)]/#2[\([1]\)]\  > \ 0.01], \ 
                  2, \ 25]], _, \ Rule]; 
data\  = \ 
          Last[data]; 
bounces\  = \ \(("Bounces"\  /. \ data)\); 
xmax\  = \ 
          Max[bounces[\([All, 1]\)]]; 
xmin\  = \ 
          Min[bounces[\([All, 1]\)]]; 
ymax\  = \ 
          Max[bounces[\([All, 2]\)]]; 
ymin\  = \ 
          Min[bounces[\([All, 2]\)]]; 
solution\  = \ 
          Function[
            Evaluate[
              Piecewise["X"\  /. \ 
                  data]]]; 
result = {"Solution" -> solution, \ "Bounces" -> bounces, \ "EndTime" -> 
              First[nw]}; 
final\  = \ "EndTime"\  /. \ 
            result; 
\(stretch[f_]\)[x_]\  := \ 
          x\  + \ f\ \((x\  - \ Mean[x])\); 
        prange\  = \ 
          Map[\(stretch[ .1]\)[{Min[#], \ Max[#]}] &, \ 
            Transpose["Bounces"\  /. \ 
                result]]; 
ramp\  = \ 
          Block[{$DisplayFunction\  = \ Identity}, 
            FilledPlot[surf[x], Evaluate[Prepend[prange[\([1]\)], \ x]], 
              opts, PlotRange -> {{xmin, \ xmax}, {ymin, \ ymax}}, 
              AspectRatio -> 
                1]]; 
bbsol\  = \ "Solution"\  /. \ 
            result; 
bbpos\  = \ 
          Function[
            Evaluate[
              Module[{bbp\  = \ bbsol[t], \ r}, 
r\  = \ 
                  Range[Length[bbp[\([1]\)]]]; 
bbp[\([1, 
                      r, 1]\)]\  = \ 
                  r; 
bbp]]]; 
bounces\  = \ "Bounces"\  /. \ result; 
frames\  = \ 
          100; 
dt\  = \ 
          final/frames; 
Table[
          Show[{ramp, 
ParametricPlot[
                bbsol[τ], {τ, \ 0, \ t}, \ 
                PlotStyle -> Blue, \ 
                DisplayFunction :> \ 
                  Identity], 
Graphics[{PointSize[
                    0.03], 
{Gray, \ 
                    Map[Point, Take[bounces, \ bbpos[t]]]}, {Black, 
                    Point[bbsol[t]]}}]}, ImageSize -> 400, 
            TextStyle -> {FontFamily -> "Arial", 
                FontWeight -> "Bold", 
                FontColor -> RGBColor[ .4,  .4,  .4]}, 
            Frame -> {True, True, False, False}, 
            FrameLabel -> {"x", "y"}, 
            AspectRatio -> Automatic], 
{t, dt, final, 
            dt}]
]\)