ODE system with discontinuous right-hand-side

Question

I have a 1st order ODE system. One of the equations is piecewise function of 2 of the dependent variables. I try to solve it in Python environment.

\begin{align}\dot x_1 &= x_2 - x_3\\ \dot x_2 &= x_1 + x_4\\ \dot x_3 &= \left\{ \begin{array}{ll} (x_2-x_{4})^{(0.5)} & x_2\ge x_{4}\\ 0, & x_2\lt x_{4} \end{array} \right.\\ \dot x_4 &= x_3 - x_1\end{align}

(Actually this is not my real ODE system. I just fabricate it to sketch my problem. Here is my real ODE system

I have read carefully the book A Primer of Scientific Programming with Python. In that book it is stated that we must write ODEs like that : \begin{align}\dot u &= f(u,t)\end{align} For scalar ODEs u and f correspond to float objects, and to arrays for system of ODEs. But didn't mention if the right hand side is piecewise function of unknown variables.

Answer 1

Most (if not all) ODE solvers have certain continuity requirements on the right-hand side. If this is not the case, an adaptive solver will usually regulate the step size smaller and smaller (and eventually fail due to hitting the minimum step size or similar). The main exception is if the integration step exactly hits the point of discontinuity.

Now, your function is only discontinuous in the derivative of the right-hand side, but this discontinuity is pretty severe (being a jump from 0 to ∞). You may try just implementing the discontinuity verbatim, but this will likely fail.

The better alternatives are:

• Smooth the discontinuity by using a very sharp sigmoid instead. This is probably also more realistic anyway.

• Use event-detection to locate the point of discontinuity and swap the right-hand side at this point.

Answer 2

Typical ODE solvers simply need to provide you a function that describes the evaluates the right hand side. That is, they will call this function with a given $x,t$ and you simply have to return what $f(x,t)$ is for these arguments. In your case, the implementation of that function would simply contain an if statement. That is a totally common thing in practice and any good ODE solver will not bat a lash over it.

The bigger issue in your case is that your right hand side is not Lipschitz. The convergence proof of most ODE solver algorithms requires that the right hand side if Lipschitz to guarantee the usual convergence order. The right hand side you have here does not satisfy this requirement -- which may or may not mean that your method of choice converges at the expected rate. But that's a question separate from what you were asking.