I am trying to numerically find an equilibrium (maximum) of a function using its differential. The following is a simplified version.
myEquilibrium[p_]:=Last[{myPreviousStep=1;NDSolve[{s'[t]==p-s[t],s[0]==myPreviousStep,WhenEvent[s[t]-myPreviousStep<10^-4||s[t]<10^-4,"StopIntegration"]},{s},{t,0,Infinity},StepMonitor:>(myPreviousStep=s[t])];myPreviousStep}]
myEquilibrium[.5] // 0.5 which is correct
The function should be constrained to positive s, which is why I included the s[t]<10^-4
requirement. However this does not seem to work.
myEquilibrium[-.5] // -0.5 but should be 10^-4
The NDSolve should also not go to 0
exactly, as the real, non-simplified version contains a 1/s[t]
in the differential. That's another reason I want the procedure to stop at s[t]=10^-4
, before s[t]=0
.
myEquilibrium[0] // 0. but should be 10^-4
Finally, I often get NDSolve::mxst: Maximum number of 10000 steps reached at the point t == 6.57563031118913721074693815431*^4952. >>
and NDSolve::ndsz: At t == 1.79769313486*^308, step size is effectively zero; singularity or stiff system suspected. >>
. What would typically solve this issue?
Thanks.