How to write down that WhenEvent?

Question

I am translating that nice animation of a ball bouncing on a hilly terrain (can be downloaded as a PDF file and as a GIF file) from Maple language into Wolfram Language. Here is my code.

surf = Sin[x] + 0.2*Cos[4 x + Sin[4 x]] - 0.2*x + 3;
Plot[surf, {x, 0, 10}, AspectRatio -> Automatic]

n = Normalize[{D[surf, x], 1}] /. x -> x[t];v = {x'[t], y'[t]};c = 0.99;
vreflect = -(1 + c)*v.n*n + v // Simplify;
deq = {y''[t] == -9.81, x''[t] == 0, x[0] == 2, y[0] == 4.5,x'[0] == 0, y'[0] == 0};
NDSolve[{Evaluate[deq], WhenEvent[Evaluate[y[t] == surf /. x -> x[t]], {y'[t] -> vreflect[[2]], 
x'[t] -> vreflect[[1]]}]}, {x[t], y[t]}, {t, 0, 10}]

However, the latest command produces an error. I need to tronslate the code of Maple events = [[y(t) = eval(surf, x = x(t)), [temp = diff(x(t), t), diff(x(t), t) = V_reflect[1], diff(y(t), t) = subs(diff(x(t), t) = temp, V_reflect[2])]]], where the derivative diff(x(t), t) is saved and used twice, into Wolfram Language.

Answer 1

The following Mathematica workflow shows that it tracks the Maple solution quite well if one decreases MaxStepFraction (other simulation parameters may also work).

Workflow

Create Bounce Function

(* Hilly surface *)
surf[x_] := Sin[x] + 1/5*Cos[4 x + Sin[4 x]] - 1/5*x + 3;
(* Calculate normal using FrenetSerretSystem *)
{tangent[x_], normal[x_]} = Last@FrenetSerretSystem[{x, surf[x]}, x];
(* Create bounce function from Maple PDF *)
vbounce[vx_, vy_, {nx_, ny_}, 
   cor_ : 0.99] = (-(1 + cor) ({vx, vy}.{nx, ny}) {nx, ny} + {vx, 
      vy}) // Simplify;

Setup and Solve ODE

(* Hill plot *)
hill = Plot[surf[x], {x, 0, 10}, AspectRatio -> Automatic];
(* Setup and solve diffeq system *)
g = -9.81;
c = 0.99;
tend = 10;
deq = {y''[t] == g, x''[t] == 0, x[0] == 2, y[0] == 4.5, x'[0] == 0, 
   y'[0] == 0};
{xfun, yfun} = 
  NDSolveValue[{deq, 
    WhenEvent[
     y[t] == surf[x[t]], {{x'[t], y'[t]} -> 
       vbounce[x'[t], y'[t], normal[x[t]], c]}]}, {x, y}, {t, 0, 
    tend}, MaxStepFraction -> 1/2000];

Import Images from Dropbox Link

(* Import gif image from Dropbox download *)
gif = Import["view.gif"];

Create Manipulate to Compare Maple and Mathematica

(* Create manipulate for comparison purposes *)
grids[min_, max_] := 
 Join[Range[Ceiling[min], Floor[max]], 
  Table[{j + .5, Dashed}, {j, Round[min], Round[max - 1], 1}]]
frames = Table[
   Rasterize@
    Show[{Graphics[{Blue, PointSize[Large], 
        Point[{{xfun[t], yfun[t]}}]}], hill}, 
     PlotLabel -> 
      Style[StringTemplate["Time = ``"][
        ToString@PaddedForm[t, {6, 4}]], 14], 
     PlotRange -> {{0, 10}, {0, 5}}, Axes -> True, 
     GridLines -> grids], {t, 0, tend, 0.067114}];
Manipulate[Column[{frames[[i]], gif[[i]]}], {i, 1, 150, 1}, 
 ControlPlacement -> Top]

Animation

(* Create Comparison Animation *)
cols = Table[Column[{frames[[i]], gif[[i]]}], {i, 1, 150}];
SetDirectory[NotebookDirectory[]]
Export["bounce.gif", cols[[1 ;; 100]], 
 "AnimationRepetitions" -> ∞]

Answer 2

The problems lies in the original code:

1. wrong tangent vector(maybe you were trying normal vector)
2. redundant or invalid Evaluate

---

I rewrote the code:

Clear@"`*"
surf[x_] := Sin[x] + 0.2*Cos[4 x + Sin[4 x]] - 0.2*x + 3;
hill = Plot[surf@x, {x, 0, 10}, AspectRatio -> Automatic]
n = Normalize@{1, surf'[x[t]]};
v = {x'[t], y'[t]};
c = 0.8;
vreflect = (1+c)*Projection[v, n] - v //Simplify;
deq = {y''[t] == -9.81, x''[t] == 0, x[0] == 2, y[0] == 4.5, x'[0] == 0, y'[0] == 0};
sol = Reap@NDSolve[
    {
        deq,
        WhenEvent[
            y[t] <= surf@x[t],
            {
                {x'[t], y'[t]} -> (Sow@{x[t],y[t]};#)
            }
        ]&@vreflect,
        WhenEvent[
            x[t] > 10 || x[t] < 0,
            tmax=t;"StopIntegration"
        ]
    }
, {x, y}
, {t, 0, 10}
, MaxStepSize -> 0.0001];
If[!ValueQ@tmax,tmax=10];
Show[
    ParametricPlot[{x@t, y@t} /. sol[[1,1]], {t, 0, tmax}, PlotStyle -> Red, PlotRange -> {{0,10},{0,4.5}}],
    hill,
    Graphics@Point@sol[[2,1]]
]

---

If you want the dissipation to be with the normal component instead of the tangent component, the formula should be:

vreflect = (1+c)*Projection[v, n] - c*v //Simplify;

But the ball will be no longer energetic. The collision will be frequent and difficult to be detect.

Answer 3

Here we provided a 3D version.

Clear["Global`*"];
reflect[vector_, 
   normal_] = -(vector - 2 (vector - Projection[vector, normal]));
f[x_, y_, z_] = z - Sin[y*Sin[x]];
g[x_, y_, z_] = x^2/4^2 + y^2/3^2 + z^2/8^2 - 1;
(*reg=ImplicitRegion[g[x,y,z]<=0,{x,y,z}];*)

reg = Ellipsoid[{0, 0, 0}, {4, 3, 8}];
surfs = ContourPlot3D[f[x, y, z] == 0, {x, y, z} \[Element] reg, 
   Mesh -> None, ContourStyle -> Darker@Cyan, 
   RegionBoundaryStyle -> Directive[GrayLevel[.1], Opacity[.1]], 
   BoxRatios -> Automatic, Axes -> False, Boxed -> False];

pt0 = {.5, .5, 6};
dir0 = {0, 0, 0};
sol = NDSolveValue[{x''[t] == 0, y''[t] == 0, 
    z''[t] == -9.8, {x[0], y[0], z[0]} == pt0, {x'[0], y'[0], z'[0]} ==
      dir0,
    WhenEvent[
     f[x[t], y[t], z[t]] == 
      0, {{x'[t], y'[t], z'[t]} -> 
       reflect[{x'[t], y'[t], 
         z'[t]}, {Derivative[1, 0, 0][f][x[t], y[t], z[t]], 
         Derivative[0, 1, 0][f][x[t], y[t], z[t]], 
         Derivative[0, 0, 1][f][x[t], y[t], z[t]]}]}], 
    WhenEvent[
     g[x[t], y[t], z[t]] == 
      0, {{x'[t], y'[t], z'[t]} -> 
       reflect[{x'[t], y'[t], 
         z'[t]}, {Derivative[1, 0, 0][g][x[t], y[t], z[t]], 
         Derivative[0, 1, 0][g][x[t], y[t], z[t]], 
         Derivative[0, 0, 1][g][x[t], y[t], z[t]]}]}]},
   {x[t], y[t], z[t]}, {t, 0, 15}, MaxStepSize -> .01];
ani = Animate[Show[surfs,
   ParametricPlot3D[sol, {t, 0, c}, PlotPoints -> 200, 
     PlotStyle -> Directive[Opacity[.5], Yellow], Mesh -> {{c}}, 
     Method -> {"BoundaryOffset" -> False}, 
     MeshStyle -> Directive[Red, Opacity[1], AbsolutePointSize[10]]] /. 
    Line -> Arrow, PlotRange -> All, Boxed -> False, 
   ImageSize -> Large], {c, $MachineEpsilon, 15}, 
  ControlPlacement -> Bottom, AnimationRate -> 1]