# Parameter which is updated in NDSolve equation

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.

• What values would you like for your constants? – bbgodfrey Oct 17 '17 at 2:24
• Try treating as a differential-algebraic equation, where [Lambda][x] is a second variable. – bbgodfrey Oct 17 '17 at 4:08
• @bbgodfrey Good point. 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. How would treating the problem as a DAE help? AFAIK this does not let me define an update rule for $\lambda$ instead of a differential relation. – Hidde Oct 17 '17 at 7:24
• In my experience, readers expect that questions be expressed in syntactically correct Mathematica code. In your case, the code immediately following "Code that works" should include expressions for Ξ and the constants, all in Mathematica form. I also recommend that you copy the code from your notebook rather than typing it into the question to avoid typographical errors. – bbgodfrey Oct 17 '17 at 13:34
• OK, I've added the working code for both situations, with and without the $\lambda$. Both situations evaluate without problems, but the second instance only uses the initial value of $\lambda$ instead of the updated value because of the StepMonitor. Thanks for your considerations. – Hidde Oct 18 '17 at 8:30

λ 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.

• Thank you for this elaborate answer. So WhenEvent[...] was the solution, where the time step can be made as small as required. Although this means that the time steps will never be exactly the same as those of the numerical integration method. Furthermore: it is known that this method (and corresponding $\Xi$) is not stable, that was exactly what I needed to visualize. – Hidde Oct 19 '17 at 10:22
• @Hidde For certain Methods It is possible to force the time-step to be a constant, in which case you could reset λ at every time-step. – bbgodfrey Oct 19 '17 at 12:32
• Yes, of course. But that is making the method conform to my requirements instead of updating the $\lambda$ according to the method's time steps. Still, your proposal works well enough for my needs. – Hidde Oct 19 '17 at 12:44