Using NDsolve with WhenEvent is not working

Question

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!

Answer 1

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

Answer 2

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