numerical solution for discountinuous differential equations

Question

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?

Thanks in advance

Answer 1

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.

Answer 2

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.