Implicit ODE solver with discontinuous derivatives

Question

I want to implement an implicit ODE solver, but don't know what to do when the differential equations (DEs) have discontinuities of the form:

1. More common type: $$\lim_{x\rightarrow 0}\frac{\sin(x)}{x} \enspace .$$

2. Less common, as it occurs in only one DE $$\min(\alpha x, \beta) \enspace .$$

I am trying to implement this myself because it is a set of massively coupled, identical DEs and we want to run them, because of our pipeline, on GPUs in Matlab. I considered adapting/modifying Matlab ode15 but I would still want exact evaluations of the Jacobian, since the system has about 10k state variables, so doing numerical approximations to the Jacobian would be expensive.

I have generally considered restarting the ODE solver after the discontinuity but I am not sure how to do that in practice: do I define an $\epsilon$ region around the discontinuity and then simply check if the solver reached it, and then restart?

For the second type of discontinuities I thought that it's the proper limit of a sequence of function of the form

$$\lim_{\gamma \rightarrow \infty}\log\left(\frac{\alpha}{1 + \exp{(-\gamma x)}}\right) \enspace ,$$

so the Jacobian would be defined piecewise as $\alpha$ for $x<0$, $0$ for $x>0$ and $\frac{\alpha}{2}$ for $x = 0$.

Answer 1

With the discussion clearing up some of the confusion in the original post, here is a summary of the discussion so far:

I want to implement an implicit ODE solver, but don't know what to do when the differential equations (DEs) have discontinuities of the form:

1. More common type: \begin{align}\lim_{x\rightarrow 0}\frac{\sin{x}}{x}\end{align}
2. Less common, as it occurs in only one DE \begin{align}\min(\alpha x,\beta)\end{align}

As Kirill points out, the first type of "discontinuity" is removable, which you can prove in this case either using the Squeeze Theorem or L'Hôpital's Rule. Without more specific information, L'Hôpital's Rule is one of the more general strategies for resolving these types of limits, so it's a good heuristic strategy to consider for other parts of your problem.

The second type of "discontinuity" is really more a point of nondifferentiability.

I am trying to implement this myself because it is a set of massively coupled, identical DEs and we want to run them, because of our pipeline, on GPUs in Matlab.

I don't know of any libraries that implement ODE solvers on GPUs in MATLAB. SimEngine is an MIT-licensed MATLAB toolbox that seems to implement an ODE solver for the GPU. It also seems to be unsupported (the company that wrote it does not appear to exist anymore, and its website no longer exists).

I don't recommend implementing this type of software yourself unless you really know what you're doing. As I've said multiple times in the comments, there are BSD-licensed libraries (such as SUNDIALS) that you might be able to modify for your purposes.

I considered adapting/modifying Matlab ode15 but I would still want exact evaluations of the Jacobian, since the system has about 10k state variables, so doing numerical approximations to the Jacobian would be expensive.

You don't necessarily need an exact Jacobian matrix, but you do need a function that supplies a good enough approximation to the Jacobian matrix using some analytical expression. (W-methods do not require an exact Jacobian matrix, and in practice, many libraries implement methods that only require sufficiently accurate Jacobians.) You are correct that evaluating the Jacobian via finite differencing would be expensive.

I have generally considered restarting the ODE solver after the discontinuity but I am not sure how to do that in practice: do I define an $\epsilon$ region around the discontinuity and then simply check if the solver reached it, and then restart?

You have two general options:

1. Ignore discontinuities, in which case your ODE solver will be first-order accurate (for a finite set of ODEs; in the infinite-dimensional case, the error will be of $O(h^{1/2})$, where $h$ is the step size).

2. If you want to use a $k$th order method, locate any discontinuities in your right-hand side or any of its derivatives up to order $k - 1$. At each of these discontinuities, restart the integrator. Typically, discontinuities are expressed in terms the roots of an "event function" that is continuously differentiable; discontinuities are then located via rootfinding. This capability is implemented in SUNDIALS, and Hairer, Nørsett, and Wanner probably discuss how to implement it in their textbook.