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$


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?

enter image description here


  • 1
    $\begingroup$ 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) $\endgroup$ Commented Oct 13, 2021 at 11:46
  • $\begingroup$ Thank you so much, i will try what you have said!!! $\endgroup$ Commented Oct 13, 2021 at 20:20

2 Answers 2


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:

![enter image description here

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

enter image description here

  • $\begingroup$ Thank you so much for the info!! $\endgroup$ Commented 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]]) == 
   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

Mathematica graphics

The initial conditions were easy to adjust by hand:

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.