4
$\begingroup$

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.

Improve this question
$\endgroup$
4
  • $\begingroup$ see this related post. you'll want to look up information on zero-crossing functions and event location for numerical ODE solvers. $\endgroup$
    – GoHokies
    Apr 30, 2018 at 13:05
  • $\begingroup$ I know that there are special solvers for this kind of systems, but what happen if you just write the code with your piecewise function using an if statement? $\endgroup$
    – nicoguaro ♦
    Apr 30, 2018 at 14:32
  • $\begingroup$ Are you asking for the basic approach to setting this up or the more specific question that the two above comments are addressing? $\endgroup$
    – Kyle Mandli
    Apr 30, 2018 at 18:18
  • $\begingroup$ @Kyle Mandli the ODE system that i gave in the question is not the actual one that I need to solve. So I want to learn basic approach and then apply it. If you look at my real ODE system , it is a little bit different. At first I will try to follow nicoguaro 's suggestion. $\endgroup$
    – kemal kaya
    Apr 30, 2018 at 18:29

2 Answers 2

Reset to default
3
$\begingroup$

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.

Improve this answer
$\endgroup$
2
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.