# Implicit ODE solver with discontinuous derivatives

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$.

• Welcome to SciComp.SE. What have you tried so far? What are the features of your DEs? Sep 7, 2015 at 18:13
• @nicoguaro, I'm a newbie to this. I've tried analytically doing the Jacobian, and obviously I run into nonsense. The DE's are non-linear and discontinuities in a finite number of points. Sep 7, 2015 at 23:09
• You are not giving us too much detail. But you can compute the Jacobian in a CAS and then export the code to your language of use. This can be done in SymPy, Maple, Mathematica and more (like Maxima to some extent). Sep 7, 2015 at 23:35
• Since it is only your ODE system that has discontinuities, you should at the very least be able to determine whether your ODE solver works correctly on its own by testing it on a different, simpler problem. There is simply not enough information in this question to give any kind of meaningful answer or a guess about what is going on here. "Run into nonsense" is not informative at all. Sep 8, 2015 at 0:04
• I generally recommend not implementing these methods yourself, except possibly on simpler problems as an exercise in learning about them. Integrating through a discontinuity is only going to reduce the order of accuracy of implicit methods. The standard solution is to integrate up to the discontinuity, stop, then reinitialize the ODE solver with the new ODE system on the other side of the discontinuity using the current solution as an initial condition. This approach is implemented in SUNDIALS. Sep 10, 2015 at 8:14

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.

• thanks for the long answer. As usual, I should have thought more before impulsively asking my question. Regarding the non-differentiability in the second problem, is it a problem at all? I'm obviously new to this, and I want some assurance that defining the Jacobian of the piecewise continuous function will not introduce errrors. Sep 11, 2015 at 12:09