1
$\begingroup$

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.

$\endgroup$
5
  • $\begingroup$ What values would you like for your constants? $\endgroup$
    – bbgodfrey
    Commented Oct 17, 2017 at 2:24
  • 1
    $\begingroup$ Try treating as a differential-algebraic equation, where [Lambda][x] is a second variable. $\endgroup$
    – bbgodfrey
    Commented Oct 17, 2017 at 4:08
  • $\begingroup$ @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. $\endgroup$
    – Hidde
    Commented Oct 17, 2017 at 7:24
  • 1
    $\begingroup$ 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. $\endgroup$
    – bbgodfrey
    Commented Oct 17, 2017 at 13:34
  • $\begingroup$ 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. $\endgroup$
    – Hidde
    Commented Oct 18, 2017 at 8:30

1 Answer 1

1
$\begingroup$

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

enter image description here

enter image description here

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]

enter image description here

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.

$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – Hidde
    Commented Oct 19, 2017 at 10:22
  • $\begingroup$ @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. $\endgroup$
    – bbgodfrey
    Commented Oct 19, 2017 at 12:32
  • $\begingroup$ 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. $\endgroup$
    – Hidde
    Commented Oct 19, 2017 at 12:44

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.