Parameter which is updated in NDSolve equation

Question

I have an NDSolve problem which is easy to solve numerically, given by

$$m y''(x)+\Xi(y(x)) = -mg, \quad y(0)=y_0, \quad y'(0)=v_0$$

where $\Xi(y(x))$ may be some non-linear function of $y(x)$.

However, now I want to add a parameter $\lambda$ which must be updated each time step. The equation becomes $$m y''(x)+\Xi(y(x),\lambda) = -mg, \quad y(0)=y_0, \quad y'(0)=v_0, \quad \lambda = \lambda_0$$ and $$\lambda^+ := \Xi(y(x), \lambda)$$ as an update rule.

I've tried multiple things:

• Making $\lambda$ a function of $x$ and making it one of the dependent variables. However, I only have an update rule and not a differential relation for $\lambda$.
• Using a StepMonitor to update $\lambda$ to a new value, each time step. This works perfectly (as far as I can see), however the differential equation given as a parameter to NDSolve is not evaluated which keeps the inital value of $\lambda$.
• The DiscreteVariables option which should allow updating a parameter. However, I see no differences in the solution or the solving process.

Any help would be greatly appreciated!

---

Code which works:

Module[{\[CapitalXi], c = 10, m = 1, g = 10, y0 = 10, v0 = 0, T = 10},
 \[CapitalXi][y_] := -c Piecewise[{{0, -y < 0}, {-y, -y >= 0}}];
 With[{sol = 
    NDSolveValue[{m y''[t] + \[CapitalXi][y[t]] == -m g, y[0] == y0, 
      y'[0] == v0}, y[t], {t, 0, T}]},
  Plot[sol, {t, 0, T}]]
 ]

Code which does not work because the equation is not evaluated again:

Module[{\[CapitalXi], p = 10, m = 1, g = 10, y0 = 10, v0 = 0, 
  T = 10, \[Lambda] = -1},
 \[CapitalXi][y_, \[Lambda]_] := 
  With[{s = 1, \[Beta] = .5}, 
   Piecewise[{{\[Lambda] 1/(s + p y), 
      y >= -((s \[Beta])/p)}, {\[Lambda] (s - p y - 2 s \[Beta]) /(
       s^2 (-1 + \[Beta])^2), y < -((s \[Beta])/p)}}]];
 With[{sol = 
    NDSolveValue[{m y''[t] + \[CapitalXi][y[t], \[Lambda]] == -m g,
      y[0] == y0,
      y'[0] == v0
      }, y[t], {t, 0, T}, 
     StepMonitor :> (\[Lambda] = \[CapitalXi][y[t], \[Lambda]])]},
  Plot[sol, {t, 0, T}]
  ]
 ]

---

For completeness, an example of $\Xi(y(x), \lambda)$ may be $$\begin{cases} \frac{\lambda }{p y(x)+s} & y(x)\geq -\frac{\beta s}{p} \\ \frac{\lambda (-p y(x)-2 \beta s+s)}{(\beta -1)^2 s^2} & y(x)<-\frac{\beta s}{p} \\ \end{cases}.$$

You may assume $p \approx 10$ (but may be increased), $s=1$, $\beta \in (0,1)$ (usually something like $0.5$. The value of $\lambda$ is negative, and $\lambda_0=-1$ is usual.

Answer 1

λ can be reset not at every time step but very often by treating it as a discrete variable and updating it periodically with WhenEvent.

Clear[λ]
p = 10; m = 1; g = 10; y0 = 10; v0 = 0; T = 10; s = 1; β = .5;
Ξ[y_, λ_] := Piecewise[{{λ 1/(s + p y), y >= -((s β)/p)}, 
    {λ (s - p y - 2 s β)/(s^2 (-1 + β)^2), y < -((s β)/p)}}]
sol = NDSolveValue[{m y''[t] + Ξ[y[t], λ[t]] == -m g, y[0] == y0, y'[0] == v0, λ[0] == -1, 
    WhenEvent[Mod[t, .1] == 0, λ[t] -> Ξ[y[t], λ[t]]]}, {y[t], λ[t]}, {t, 0, T}, 
    DiscreteVariables -> {λ}, MaxStepSize -> 0.01];
Plot[First@sol, {t, 0, T}, AxesLabel -> {t, y}, 
    LabelStyle -> Directive[Black, Bold, Medium], ImageSize -> Large]
Plot[Last@sol, {t, 0, T}, PlotRange -> {-.1, .1}, AxesLabel -> {t, λ}, 
    LabelStyle -> Directive[Black, Bold, Medium], ImageSize -> Large]

To understand the behavior of λ[t], it may be helpful to plot Ξ[y, λ], which is the reset value of λ.

Plot3D[Ξ[y, λ], {y, -1, 1}, {λ, -1, 1}, AxesLabel -> {y, λ}, PlotRange -> All, 
    LabelStyle -> Directive[Black, Bold, Medium],ImageSize -> Large]

When y > -(s β)/p, here equal to -0.05, λ shrinks in size by a few orders of magnitude for each reset. For instance, the first reset of the computation,

Ξ[y0, -1]
(* -(1/101) *)

abruptly reduces λ from -1 to -1/101. On the other hand, when y < -(s β)/p each reset increases λ by a few orders of magnitude in size. Hence, the solution appears to be dependent on the size of the initial time steps. Perhaps, this problem is not well posed.