# numerical solution for discountinuous differential equations

Let a system of differential equations has the following constraint:
$$\frac{d^2x}{dt^2}=f_x(t, x,\frac{dx}{dt})$$ $$\frac{d^2y}{dt^2}=f_y(t, y,\frac{dy}{dt})$$ if $$\frac{dx}{dt} > \frac{dy}{dt}$$ and $$\frac{d^2x}{dt^2}=g_x(t, x,\frac{dx}{dt})$$ $$\frac{d^2y}{dt^2}=g_y(t, y,\frac{dy}{dt})$$ otherwise.

Obviously, There is a hard nonlinearity in this system of equations. Actually, when two rigid bodies collide, this kind of change in the governing equations may happen. The worst is that the condition is defined based on the "speeds" here, not on the "positions'.
Furthermore, in this case, the discontinuity is not in the right-hand side. But is in the dynamic of system itself.
The implicit numerical algorithms (for example: implicit Runge Kutta's) seem not to help here; Maybe, because they work with the end-point data of the time step and due to the conditional behavior of this system of equations, we cannot rely on the end-step data.

The question is: Is it possible to solve this system by the MATLAB ODE functions? If so, which algorithm is more reliable? Is it possible to feed the system of equations along with its condition entirely to the algorithm (as a odefunction) without any intervention?

I am pretty sure you are looking for the Matlab documentation about ODE Event Location.

The fact that the condition is based on speed shouldn't be a problem. You need to rewrite your 2nd order equations into a system of first order equations: \begin{equation} \begin{aligned} \dot{x}_1 &= v_1\\ \dot{v}_1 &= f_{x} (t, x_1, v_1)\\ \dot{y}_2 &= v_2\\ \dot{v}_2 &= f_{y} (t, x_2, v_2) \end{aligned} \end{equation} with added information about speed being derivative of position. Your discontinuity than turns into: \begin{equation} v_2 - v_1 < 0 \end{equation} This equation you can use to stop the integration and change the equations using the logic similar to bouncing ball and restricted three body problem examples in the above documentation. The documentation tells you to type

ballode

in the Matlab command line in order to run the example. This just shows you the output of the integration. In order to see how this is implemented I suggest you check the code of the example. If you type

edit ballode

in the Matlab command line you will get the code in the text editor. For the bouncing ball it iterates in a for loop where each iteration of the loop represents an integration of the system until next event. At the beginning or the of each iteration you can switch system to represent the change of the behavior.

• You mean that MATLAB automatically detects the point where this inequality satisfied, stop integrating, and then continues the process as a new initial value problem from this point? What if this "event" happens several times during the simulation? Sep 14, 2018 at 10:01
• Another question: can MATLAB switch between the two sets of differential equations? The problem here is that we have to replace the differential equations, not only the "right-hand side" of them. Sep 14, 2018 at 10:16
• I edited the answer to add the command which lets you see the code for the bouncing ball example. You can probably do a copy/paste and adapt it to your system. Admittedly, your system is more complicated, but I think that modifying the given example can save you a lot of time, rather trying to implement everything by yourself. Sep 14, 2018 at 15:33
• I wrote a odefunction for this problem in the form of "conditional function", exactly in the form that I described the problem above (except that I represented it in the state-space form, as you mentioned in your note) The Matlab ode codes were successful to evaluate the equation a bit. But the answers were not reliable for all setups of the problem. I think some manipulations are necessary. I think it is better to intervene the solution process and do not rely on the Matlab functions entirely, as Yves advised. Thanks for the good help Sep 15, 2018 at 20:04

I agree with @slavenglumac's answer, there's nothing really special here, except that the functions $f/g$ have an extra dependency on the sign of $\dot x-\dot y$. They can be merged into one.

If you want a better handling of the discontinuity, you can integrate until you observe a change of sign, and estimate the zero by interpolation. Then restart from this point.