# How to write down that WhenEvent?

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.

• The downvoter: Can you ground your vote? TIA. Sep 1, 2020 at 13:06
• I'm not the downvoter, but currently the question doesn't look like a self-contained one i.e. it seems that answerers have to read an external PDF file (involving Maple code! ) to understand the whole question. Explaining the system with traditional math notation first and then showing the unworking code may be better. Sep 4, 2020 at 2:25

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 = "][
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" -> ∞]


• +1. Many thanks from me to you for your work. You are a Mathematica master! First, do you mean ODE instead of PDE? Second, gif = Import["view.gif"] does not work for me (perhaps, some path should be added). Third, can you extend both the animations to 20 sec? Fourth, a good code is a commented code (Comments are useful for the author as well for the reader.). In order to approve your answer, can you kindly respond my requests? Sep 2, 2020 at 16:11
• gif = Import[ https : // www.dropbox.com/s/q9mnz46tyf2az77/view.gif?dl = 0, "view.gif"] does not work too. Sep 2, 2020 at 16:48
• If I extend the time in the Maple application to t = 0 .. 150, then the ball falls. Also ( if I am not mistaken) the result of ...{t, 0, tend, 0.5}]; Manipulate[Column[{frames[[i]]}], {i, 1, 150, 1}, ControlPlacement -> Top] with tend = 150; differes from the Maple's one in this case Sep 2, 2020 at 17:52
• @user64494 I had to download the gif file manually from the Dropbox and I put it in the saved notebook directory. You probably could transfer the file to your own personal Dropbox location and import using Insert>FilePath. The way I synched the two simulation was I counted the number of frames in the gif animation and I set the increment to timestamp I observed in gif[[2]]. You will have to do the same to sync the timestamps. Also note, the gif animations are limited to 2MB on this site, which limits how many frames you can show. Sep 2, 2020 at 19:02
• @user64494 I think that you are substantially moving the goal post from your initial question from getting a particular simulation to work to one that is a general robust solution to a large parameter space. I agree with @xzczd's comment about MaxStepSize since you are interested in increasing tend 15x. You should also look at the Possible Issues>Arbitrarily close events, which happens when the particle slows down to near zero velocity. To get $tend=150$ to work, try MaxStepSize -> 1/200, WorkingPrecision -> 50, MaxSteps -> Infinity Sep 3, 2020 at 19:40

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.

• I wonder the upvoter, though here mathematica.stackexchange.com/questions/223847/… I see 7 votes for the wrong answer. Sep 1, 2020 at 11:43
• Here is my comment, deleted by an unknown to me person. Thank you for your work. A good code is a commented code. There are comments in the Maple application, but there aren't comments to your code. First, this is not an animation. Second, your result differs from the result in the application. Third, your "wrong tangent vector" is ungrounded. Sep 2, 2020 at 4:44
• Put c=0.8 in your code and explain the result. TIA. Sep 3, 2020 at 0:19
• @user64494 Both of you have defined the normal vector in wrong way… Just plot it: origin = {x, surf}; nvec = Normalize[{D[surf, x], 1}]; arrow[x_] = {origin, origin + nvec}; Manipulate[plot~Show~Graphics[Arrow@arrow@x], {x, 0, 10}]. The correct normal vector is nvec = Normalize[{-D[surf, x], 1}];. Sep 4, 2020 at 1:55
• @user64494 If you're confused about the {x'[t], y'[t]} -> vreflect, this is discussed in Possible Issues section of document of WhenEvent, starting from the line With sequential event actions, the variables are modified in turn…. Sep 5, 2020 at 5:40

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]