# Using NDsolve with WhenEvent is not working

I am trying to use WhenEvent in this nonlinear equation system:

$$\ddot r-\dot\phi^2r+9.8(\sin\phi-0.7\cos\phi)=0$$

$$1.2(1+20(r^2+2^2+4r\sin\phi)^2)\ddot\phi+\dot\phi+\phi=\sin4t$$

I want that When $$r(t)>1 \Rightarrow r(t)=1$$ and When $$r(t)<-1 \Rightarrow r(t)=-1$$, so I was using the code below:

sol = NDSolve[{r''[t] - phi'[t]^2 r[t] +
9.8 (Sin[phi[t]] - 0.7*Cos[phi[t]]) == 0,
1.2 (1 + 20*(r[t]^2 + 2^2 + 4 r[t] Sin[phi[t]])^2) phi''[t] +
phi'[t] + phi[t] == Sin[4 t],
phi[0] == phi'[0] == r[0] == r'[0] == 0,
WhenEvent[r[t] > 1, r[t] -> 1],
WhenEvent[r[t] < -1, r[t] -> -1]}, {r, phi}, {t, 0, 20}]

Plot[r[t] /. sol, {t, 0, 20}]


But when I run the code and plot $$r$$, I get values bigger than 1... What am I doing wrong?

Thanks!

• You can "fix" this by setting WhenEvent[r[t]>1,r[t]->0.9998]. This is necessary due to how WhenEvent triggers. However, DO NOT DO THIS. Since your DE is of the second order, it will still tend towards the one anyway since its "velocity" keeps unchanged. Thus you'd just create a very fine sawtooth pattern. So you'd have to set the first derivative as well to get something stable. Something like WhenEvent[r[t]>1,{r[t]->0.99999,r'[t]->0}] or even the reflection. (Add PlotRange->All to your plot here) Oct 13, 2021 at 11:46
• Thank you so much, i will try what you have said!!! Oct 13, 2021 at 20:20

We have a misunderstanding here. r[t]>1 means: the time when r[t] crosses 1 from below and NOT when r[t] is larger than 1. Look it up in the help.

To achieve what you want, that r[t] stays constant after hitting 1, you may define a piecewise function:

sol = NDSolve[{1.2 (562500 +
1.2*(r[t]^2 + 40^2 + 2*r[t]*40*Sin[phi[t]])) phi''[t] +
1030400*phi'[t] + 0.41*phi[t] == 2*Sin [4 t],
r''[t] - 0.7*9.8*Cos[phi[t]] + 9.8 Sin[phi[t]] == 0,
phi[0] == r[0] == phi'[0] == r'[0] == 0,
WhenEvent[r[t] < -1, r[t] -> -1],
WhenEvent[r[t] > 1, {r[t] -> 1, tend = t, "StopIntegration"}]
}, {r[t], phi[t]}, {t, 0, 20}]
fun[t_] =
Piecewise[{{Evaluate[r[t] /. sol[[1, 1]]], t < tend}, {1, True}}]


This gives the following interpolating function:

Plot[fun[t], {t, 0, 20}]


• Thank you so much for the info!! Oct 13, 2021 at 20:20

This sort of problem (autonomous) can be handled by constructing the desired vector field instead of using WhenEvent, since the discontinuity (event) depends on the state vector $$(r, r', \phi, \phi')$$ is independently of what time $$t$$ is. One can think of it as where versus when. It is possible to handle the problem with WhenEvent: wherever the state of the system, it happens some time or another. The issue with the OP's WhenEvent is that it does not change the vector field, and as Julien Kluge pointed out, the solution continues to increase from r[t] == 1 but does not cross it from below, as Daniel Huber pointed out.

In nonautonomous systems or events, one can use DiscreteVariables to turn on/off components of systems as shown in the documentation. In the OP's problem, we can use Piecewise to zero out the r, r' components of the vector field. NDSolve will handle the discontinuity for us. We need to convert the r components to a first-order system, so that we can annihilate both r'[t] and r''[t].

It's probably easiest to write the Piecewise form by hand, but since the OP's system wasn't in r''[t] ==... form and my eyes are unreliable, doing it by code was my preference. (It also helps that I knew how to do it by code in a couple of ways.)

ode = { (* OP's ode system *)
r''[t] - phi'[t]^2 r[t] + 9.8 (Sin[phi[t]] - 0.7*Cos[phi[t]]) ==
0,
1.2 (1 + 20*(r[t]^2 + 2^2 + 4 r[t] Sin[phi[t]])^2) phi''[t] +
phi'[t] + phi[t] == Sin[4 t]};

ode1 = Solve[
ode~Join~NestList[D[#, t] &, s[t] == r'[t], 1],
{phi''[t], s'[t], r'[t]}, {r''[t]}] /.  (* eliminates r''[t] *)
HoldPattern[d : s'[t] | r'[t] -> v_] :>
d -> Piecewise[{{v, -1 < r[t] < 1}}] /. Rule -> Equal;
ode1 // Transpose // Grid


The initial conditions were easy to adjust by hand:

sol = NDSolve[{
ode1,
phi[0] == phi'[0] == r[0] == s[0] == 0}, {r, phi}, {t, 0, 20}];

Plot[r[t] /. sol, {t, 0, 20}]