Differential Equations and Events

When in the Course of human events ...

Rob Knapp

Wolfram Research, Inc.

Initialization

Off[General::spell];
Off[General::spell1];

Needs["GUIKit`"]

Needs["NumericalMath`EquationTrekker`"]

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]

Needs["DifferentialEquations`NDSolveProblems`"]

abcsystem = 
        GetNDSolveProblem["ArnoldBeltramiChildress"];
abceqns  =  abcsystem@"System"[];
abcvars  =  abcsystem@"DependentVariables"[];

MaxTime[ifun_InterpolatingFunction]  :=  
      Max[First[InterpolatingFunctionDomain[ifun]]]

\(\(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};

\(BouncingBall[surf_, \ x0_, \ xp0_, \ 
        opts___]\  := \ 
Module[{flops, \ ν, \ k, \ 
          data, bounces, \ solution, \ nw}, 
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)\); 
solution\  = \ 
          Function[
            Evaluate[
              Piecewise["X"\  /. \ 
                  data]]]; 
{"Solution" -> 
            solution, \ "Bounces" -> 
            bounces, \ "EndTime" -> 
            First[nw]}
]\)

e·vent  
1. a happening or occurrence. esp. when important
2. a result; consequence; outcome
3. a particular contest or item in a program
-- from Webster's New World Dictionary

e·vent  
1. a happening or occurrence. esp. when important
2. a result; consequence; outcome
3. a particular contest or item in a program
-- from Webster's New World Dictionary

An event in a differential equation occurs at values of the (timelike) independent variable where some condition depending on the problem variables is satisfied.

e·vent  
1. a happening or occurrence. esp. when important
2. a result; consequence; outcome
3. a particular contest or item in a program
-- from Webster's New World Dictionary

An event in a differential equation occurs at values of the (timelike) independent variable where some condition depending on the problem variables is satisfied.

For the ODE X' (t)  ==  f(t,  X(t)) an event can be described as occurring when g(t, X(t))  =  0 for some real-valued function g, or when B(t,  X(t)) changes from True to False or vice-versa for a Boolean valued function B.

A simple example of using event detection is to find the period of a regular oscillation.  Here I will show you how to do it for the simple pendulum equation using NDSolve.

A simple example of using event detection is to find the period of a regular oscillation.  Here I will show you how to do it for the simple pendulum equation using NDSolve.

In NDSolve, the EventLocator method is a built-in controller method which detects events and takes action when they are detected.

This integrates the pendulum equation until the event condition x' (t)  =  g(t,  x(t),  x' (t))  =  
          0 is satisfied after the initial condition.  The default action to take when an event is detected is to stop the integration.

hsol  =  
      First[NDSolve[{x''[t]  +  Sin[x[t]]  ==  0,  
            x[0]  ==  3,  x'[0]  ==  0},  
          x,  {t, 0, ∞},  
          Method -> {"EventLocator",  "Event" -> x'[
                  t]}]]

A convenient way to get the final time out of the InterpolatingFunction solution returned is to use the InterpolatingFunctionAnatomy utility package.

This loads the InterpolatingFunctionAnatomy utility package.

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]

This defines a function which gets the largest time from an InterpolatingFunction. (Assuming that "time" is the first variable in the multidimensional case.)

MaxTime[ifun_InterpolatingFunction]  :=  
      Max[First[InterpolatingFunctionDomain[ifun]]]

This gets the final time and plots the solution in the phase plane up to that time.

end  =  MaxTime[x  /.  hsol];
ParametricPlot[{x[t], x'[t]}  /.  hsol, {t, 0, 
          end}];

The integration has only gone over the first part of the orbit.  In this case, by symmetry, this is half of the period.  A way to set the event at the full period only is to specify the direction of zero crossing you are interested in.

This integrates the pendulum equation for a full period.  The event is the condition x' (t)  =  g(t,  x(t),  x' (t))  =  
          0 combined with the additional condition that sign(((∂g(t,x(t),x'(t)))/(∂t))) = 1.

psol  =  
      First[NDSolve[{x''[t]  +  Sin[x[t]]  ==  0,  
            x[0]  ==  3,  x'[0]  ==  0},  
          x,  {t, 0, ∞},  
          Method -> {"EventLocator",  "Event" -> x'[
                  t],  "Direction" -> -1}]]

This plots the solution over the time interval found.

ParametricPlot[{x[t], x'[t]}  /.  psol, {t, 0, 
          MaxTime[x  /.  psol]},  AspectRatio -> 1];

Using this , it is quite simple to write a function which computes the period for a given starting angle (Subscript[x, 0]).

This defines a function which uses event location to find the period of a pendulum started from an angle x0.

\(PendulumPeriod[x0_]\  := \ 
      Module[{t, \ x, \ end, \ psol}, 
psol\  = \ 
          First[NDSolve[{\(x''\)[t]\  + \ Sin[x[t]]\  == \ 0, \ 
                x[0]\  == \ x0, \ \(x'\)[0]\  == \ 0}, \ 
              x, \ {t, 0, ∞}, \ 
              Method -> {"EventLocator", \ "Event" -> \(x'\)[t], \ "Direction" -> \(-1\)}]]; 
MaxTime[
          x\  /. \ psol]]\)

This makes a plot of the period of the pendulum as a function of starting angle.

Plot[PendulumPeriod[(Subscript[x, 0])], {(Subscript[x, 0]),  0,  π}];

Some events are expressible only by Boolean conditions as opposed to in terms of a continuously varying event function.

An example of this would be a button for stopping an integration in an interactive tool.

This opens up a small window with a button, which when pressed changes the value of the variable stop to True from its initialized value of False.

NDSolve`stop  =  False;
\(\(GUIRun[
        Widget["Panel", \ {Widget["Button", {
"label" -> "Stop Already", 
BindEvent["action", 
Script[
                    NDSolve`stop\  = \ True]]}]}]];\)\)

This integrates the Lorenz equations with an extrapolation method until the button is pressed.  Note that the event is specified with RuleDelayed so that the symbol will be evaluated each time the event is evaluated during the integration.

NDSolve`stop  =  False;
sol=NDSolve[{
x'[t]==-3 (x[t] - y[t]),x[0]  ==  0,
y'[t]==-x[t] z[t] + 26.5 x[t] - y[t],y[0]  ==  1,
z'[t]==x[t] y[t] - z[t],z[0] == 0},
{x, y, z},{t, ∞},MaxSteps -> ∞,
Method->{"EventLocator","Event":>NDSolve`stop}]

This makes a plot of part of the solution computed for the Lorenz equations.

times  =  
      Take[InterpolatingFunctionCoordinates[
            x  /.  First[sol]][[1]],  10000]; 
    Show[Graphics3D[
        Line[Transpose[{x[times],  y[times],  z[times]}  /.  
              First[sol]]]]];

The problem with this example is that NDSolve is just too darn fast.  By the time I'm able to stop the integration, the solution is too big to plot very quickly.

For the Lorenz equations, one way to look at the long term behavior along the (strange) attractor is to compute a sort of symbolic dynamics.  One possibility along this line is to record a symbol for each loop the trajectory makes: symbol 0 on the loop above the third quadrant and symbol 1 on the loop above the first quadrant.

Doing this requires that we detect each loop the trajectory makes.  An easy way to do this is to set an event at the top of a loop: z'(t) == 0 (direction -1)  Looking at the plot of the trajectory, you can see that you can detect which loop you are on by checking the signs of x(t) and y(t).  You don't want the integration to stop when these events are detected, however, you want to just save the symbol and continue on.

Reap and Sow provide a very good way to collect data for which the length is a-priori unknown in Mathematica.

Reap[Sow[1];  Sow[b];  Table[Sow[i^2], {i, 4}]]

It is possible to combine multiple events with the EventLocator method.  

NDSolve[..., Method->{"EventLocator", "Event"->{ (Subscript[e, 1]),  (Subscript[e, 2]),   ...  ,  (Subscript[e, n])}, 
									  "Direction"->{ (Subscript[d, 1]),  (Subscript[d, 2]),   ...  ,  (Subscript[d, n])}, 
									  "EventAction"->{ (Subscript[a, 1]),  (Subscript[a, 2]),   ...  ,  (Subscript[a, n])}}]

detects events (Subscript[e, 1]),  (Subscript[e, 2]),   ...  ,  (Subscript[e, n]) with corresponding direction restrictions (Subscript[d, 1]),  (Subscript[d, 2]),   ...  ,  (Subscript[d, n]) and actions (Subscript[a, 1]),  (Subscript[a, 2]),   ...  ,  (Subscript[a, n]).   If several of the events occur in a single step, the actions are always evaluated in integration order.  If only one direction restriction or action is given, it is used for all the events.

Reap and Sow provide a very good way to collect data for which the length is a-priori unknown in Mathematica.

Reap[Sow[1];  Sow[b];  Table[Sow[i^2], {i, 4}]]

Reap and Sow are well suited for collecting the results of event actions since you typically don't know how many times a particular event may be encountered.  The way to do this is to use Sow in the event action and Reap wrapping the use of NDSolve.

This computes the symbols described above for the Lorenz equations until the integration is stopped by pushing the button.

NDSolve`stop  =  False;
reaps = Reap[NDSolve[{
x'[t]==-3 (x[t] - y[t]),x[0]  ==  0,
y'[t]==-x[t] z[t] + 26.5 x[t] - y[t],y[0]  ==  1,
z'[t]==x[t] y[t] - z[t],z[0] == 0},
{},{t, ∞},MaxSteps -> ∞,
\(Method -> {"EventLocator", 
"Event" :> {x[t]\ y[t]\  - \ z[t], 
                    NDSolve`stop}, \ 
"Direction" -> {\(-1\), \ 
                    All}, \ 
"EventAction" :> {Sow[\((1\  + \ Sign[x[t] + y[t]])\)/2], Throw[Null, \ "StopIntegration"]}}\)]]

This forms the symbols into a nearly square matrix and makes a simple visualization of the 1's and shows it next to a similar plot of a random matrix.

symb  =  reaps[[2, 1]];
n  =  Length[symb];
n  =  Floor[Sqrt[n]];
Show[
        GraphicsArray[
          Block[{$DisplayFunction  =  
                Identity},  {{ArrayPlot[Partition[symb,  n]],  
                ArrayPlot[
                  Table[Random[Integer], {n}, {n}]]}}]]];

"The idea of reducing the study of continuous time systems (flows) to the study of an associated discrete time system (map) is due to Poincaré who first utilized it in his studies of the three body problem in celestial mechanics". ... "Nowadays virtually any discrete time system that is associated with an ODE is referred to as a Poincaré map."  [Wiggins, Inroduction to Applied Nonlinear Dynamical Systems and Chaos, 1990]

There are a numerous advantages to using a Poincaré map  to study a dynamical system. Among the most important are dimensional reduction and conceptual clarity: often results can be shown for the Poincaré map which cannot so easily be proven for the continuous system.

Often, a Poincaré map is the return map for a section (or cutting plane) across a periodic orbit.  Appropriate choice of the cutting plane will depend on the dynamical system and properties being studied.  In any case, the intersection of an orbit with a cutting plane can be found numerically by expression the cutting plane as an event.  Sometimes such Poincaré maps are referred to as Poincaré sections.

One interesting example of a Poincaré map is for the Arnold-Beltrami-Childress (ABC) flow, a widely studied model for
volume-preserving three-dimensional laminar flows used to model, among other things, solar dynamos.

This loads a collection of example problems for NDSolve.

Needs["DifferentialEquations`NDSolveProblems`"]

This gets the equations and variables associated with the ABC flow ODE system.

abcsystem = 
        GetNDSolveProblem["ArnoldBeltramiChildress"];
abceqns  =  abcsystem@"System"[];
abcvars  =  abcsystem@"DependentVariables"[];

This loads EquationTrekker which allows interactive investigation of differential equations.

Needs["NumericalMath`EquationTrekker`"]

This begins an interactive investigation of the ABC flow system showing the Poincaré map generated by the cutting plane (Subscript[Y, 2])  =  0.

\(EquationTrekker[{abceqns, \ Y\_2[0] == \ 0}, 
      abcvars, \ {T, 0, 1000}, 
      TrekGenerator -> {PoincareSection, \ "SectionCondition" -> Y\_2[T], \ "SectionVariables" -> {Y\_1[T], 
              Y\_3[T]}, \ "FilterFunction" -> 
            Function[Mod[#, \ 2  π]], 
          MaxSteps -> 
            Infinity}, 
PlotRange -> {{0, 
            2\ π}, {0, 2  π}}]\)

The Poincaré map is computed in EquationTrekker using NDSolve.  
The actual NDSolve command given uses EventLocator.

This makes a plot for a particular initial condition.  The NDSolve command used is shown in red text with the EventLocator method specification shown in blue text.

data = Reap[NDSolve[{abceqns,  (Subscript[Y, 2])[0]  ==  0,  
                  (Subscript[Y, 1])[0]  ==  2,  (Subscript[Y, 3])[0]  ==  2},{},{T, 0, 1000}, Method->\({"EventLocator", \ 
"Event" -> Y\_2[T], \ "EventAction" :> 
                      Sow[Mod[{Y\_1[T], \ Y\_3[T]}, \ 
                          2\ π]], \ "EventLocationMethod" -> "LinearInterpolation"}\), 
MaxSteps -> ∞]];
ListPlot[data[[2, 1]]];

The EventLocator method works by monitoring each step taken by the method it controls for a sign (or parity) change in any of the event functions.  When a change is detected, then the actual location of the event within the step is refined by an event location method.

The default event location method is to use FindRoot, which uses Brent's method for bracketed roots to tolerances consistent with the local truncation error.   When dense output is available from the method, it is used to compute the solution at the intermediate values when FindRoot needs to evaluate the event function.  If not, for one-step methods, the intermediate solution can  be found by taking a step of the method of appropriate size.

Sometimes, refinement of the event location to local tolerances is not necessary or desirable.  For example, when the stop button is pressed, it is actually meaningless to try to refine further.  Also, for Poincaré maps, while it is important to compute the solutions accurately to avoid accumulating errors, it is typically not possible to graphically see differences on the order of the default tolerances.

EventLocator has an "EventLocationMethod" option which allows you to specify how the event location is refined.

In cases where the event begins a substantial change in dynamics, it may be necessary to use NDSolve multiple times.  An example where it is best to use NDSolve several times is the case of the classic bouncing ball problem.

The dynamics of the ball are governed by the simple equations for a body in motion: m a  =  F.  
The forces on our idealized ball with mass 1 are gravity and wind resistance (taken to be on the order of the squared velocity.  This can be described by 

X''(t) = -ν ‖X'(t)‖ X'(t) - (Overscript[g,  ->]) 

for a position X with gravity vector (Overscript[g,  ->])  and wind drag coefficient ν.

When the ball hits the surface, we assume the bounce is a reflection but with loss of momentum determined by the coefficient of restitution k, such that the momentum immediately after the bounce is k times the momentum before the bounce (0 < k< 1).

Here is an example in 2-dimensions.

This computes the bouncing ball solution for a ball dropped onto a wavy ramp.

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

This shows a sequence of plots which can be used to animate the solution dynamics.

end  =  "EndTime"  /.  bb;
stretch[f_][x_]  :=  x  +  f (x  -  Mean[x]); 
    prange  =  
      Map[stretch[ .1][{Min[#],  Max[#]}] &,  
        Transpose["Bounces"  /.  bb]];
ramp  =  
        Block[{$DisplayFunction  =  Identity}, 
          Plot[wavyramp[x], Evaluate[Prepend[prange[[1]],  x]], 
            PlotRange -> prange, 
            AspectRatio -> 1]];
bbsol  =  "Solution"  /.  bb;
bbpos  =  
        Function[
          Evaluate[
            Module[{bbp  =  bbsol[t],  r}, 
r  =  
                Range[Length[bbp[[1]]]]; 
bbp[[1, r, 
                    1]]  =  
                r; 
bbp]]];
bounces  =  "Bounces"  /.  bb;
frames  =  100;
dt  =  end/frames;
Do[Show[{ramp, 
ParametricPlot[
            bbsol[τ], {τ,  0,  t},  
            PlotStyle -> RGBColor[1, 0, 0],  
            DisplayFunction :>  
              Identity], 
Graphics[{PointSize[
                0.04], 
{RGBColor[0, 1, 0],  
                Map[Point, Take[bounces,  bbpos[t]]]}, {RGBColor[0, 0, 1], 
                Point[bbsol[t]]}}]}], 
{t, dt, end, 
        dt}]

Here is an example in 3-dimensions

This defines a functions which represents a rippled surface and shows a surface plot of it.

ripples[x_,  y_]  =  
        Module[{r  =  x^2  +  y^2},  
          r  +  Cos[π r]/( .5  +  
                  r)];
surface  =  
        Plot3D[ripples[x, y], {x, -2.5, 2.5}, {y, -2.5, 2.5},  
          PlotPoints -> 50,  Mesh -> False];

This computes the path of a ball tossed above the surface.

bb  =  BouncingBall[ripples, {-1, 1, 8}, {2, 2, 0}]

This shows a sequence of plots which can be used to animate the solution dynamics.

end  =  "EndTime"  /.  bb;
stretch[f_][x_]  :=  
        x  +  f (x  -  Mean[x]);
bbsol  =  "Solution"  /.  bb;
bbpos  =  
        Function[
          Evaluate[
            Module[{bbp  =  bbsol[t],  r}, 
r  =  
                Range[Length[bbp[[1]]]]; 
bbp[[1, r, 
                    1]]  =  
                r; 
bbp]]];
bounces  =  "Bounces"  /.  bb;
frames  =  50;
dt  =  end/frames;
curve   =  Take[bounces,  1];
Do[curve  =  
        AppendTo[curve,  
          bbsol[t]]; 
Show[{surface, 
Graphics3D[{{RGBColor[1, 0, 0], 
                Line[curve]}, {PointSize[
                  0.03], 
{RGBColor[0, 1, 0],  
                  Map[Point[#  +  {0, 0,  .025}] &, 
                    Take[bounces,  bbpos[t]]]}, {RGBColor[0, 0, 1], 
                  Point[bbsol[t]]}}}]}], 
{t, dt, end, 
        dt}]

Annotated code which computes the solution for the bouncing ball problem is shown below.  The actual code which has been entered is in the initialization for this notebook.

The function OneBounce computes the solution from an initial condition (position X(t) and velocity X'(t)) up to the point of a bounce, returning the time and position of the bounce (on the surface surf) and the velocity immediately after the bounce reflection.  The function is designed to work with points in 2- or 3-dimensions equally well.

OneBounce[surf_, k_, ν_][{t0_,x0_, xp0_}] := 
  Module[{sol, t1, x1, x1p, g},

The gravity vector (Overscript[g,  ->]) .  This is done as a function with only numerical evaluation to avoid threading prior to numerical evaluation.

    g[t_?NumberQ] =Module[{gv = 0. x0}, gv[[-1]] = -9.8; gv];

Computes the solution X''(t) = -ν ‖X'(t)‖ X'(t) - (Overscript[g,  ->])  up until the event that the trajectory intersects the surface.  MaxStepSize is restricted so that for quickly varying surfaces, the trajectory does not pass through the surface and back through again with the event being missed.

    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]];

Set t1, x1 and xp1 to the time, position, and velocity right before the bounce reflection. Then reset xp1 to the velocity just after the bound reflection.

    t1 = InterpolatingFunctionDomain[x /. sol][[1,-1]];
    {x1, xp1} = {x[t1], x'[t1]} /. sol;
    xp1= Reflection[surf, k][{x1, xp1}];

Sow data for the solution and return.  The "X" data will be used to form a Piecewise function.

    Sow[{x[#] /.sol,t0 ≤ # ≤  t1}, "X"];
    Sow[x1, "Bounces"];
    {t1, x1,xp1}]

The function Reflection computes the reflection.    The reflection is modelled on reflection off of the tangent plane.  The velocity parallel to the plane remains the same.  Only the velocity normal to the plane is reflected.  This is easily computed with the projection against the normal.  

Reflection[surf_, k_][{x_, xp_}] := 
  Module[{zeroes, gsurf, normal, xpnew, n = Length[x] - 1},
  
The normal is computed by the gradient of z - f = c, which is (-∇f, 1).  Since the surface is given as a function without arguments (e.g. a pure function), we compute ∇f using Derivative.  We check to make sure it was computed numerically.  The normal is normalized to the unit normal vector to make the projection computation easier.

    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];
   
We don't allow glancing hits on the surface.  This could be done, but would require different dynamics since the ball would typically roll along the surface.

   If[ normal.xp == 0,
      Print["No reflection"];
      Throw[$Failed]];
   
Compute the new velocity by subtracting twice the component normal to the surface (one to cancel it and one to negate).  The result is multiplied by k for the momentum reduction.

    xpnew = k (xp - 2 normal normal.xp)
    ]

The main interface function is BouncingBall, which handles a sequence of bounces.

It has options which allow controlling of the wind drag and restitution coefficients.

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

As arguments, it takes a function (without arguments) specifying the surface, starting point (Subscript[X, 0]) and velocity (Subscript[X, 0])'.

BouncingBall[surf_, x0_, xp0_, opts___] := 
  Module[{flops, ν, k, data,bounces, solution, nw},

Check the values of the drag and restitution coefficients.

    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]];

Make sure the starting position is above the surface!

    If[Apply[surf, Most[x0]] ≥ Last[x0],
      Print["Start above the ramp"];
      Return[$Failed]];
 
Iterate bounces until one fails or until the successive time between bounces changes by less than 1% (from Shampine).  Reap gets the data saved in the OneBounce function with Sow.

   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);
  
This forms a Piecewise function that will correctly give the solution between bounces.

   solution = Function[Evaluate[Piecewise["X" /. data]];
    {"Solution"->solution, "Bounces"->bounces, "EndTime"->First[nw]}
   ]

\(BouncingBall[surf_, \ x0_, \ xp0_, \ 
        opts___]\  := \ 
Module[{flops, \ ν, \ k, \ 
          data, bounces, \ solution, \ nw}, 
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)\); 
solution\  = \ 
          Piecewise["X"\  /. \ 
              data]; 
{"Solution" -> 
            solution, \ "Bounces" -> 
            bounces, \ "EndTime" -> 
            First[nw]}
]\)

EventLocator is a controller method for NDSolve which provides an easy way to stop integration or take some other action when an event is detected.