Solving ODE with "Jumpy" Coefficients

I'm numerically solving a linear coupled ODE of the form $$y^{\prime}(t) = \hat{M}(t)y(t)=\left[\begin{array}{cc}0& A(t)\\ B(t)& 0\end{array}\right]y(t),$$ and the difficulty I'm running into is that $$A(t)$$ and $$B(t)$$ don't change for long periods of time and suddenly change rapidly. Here is a plot of one example I'm dealing with:

Blue line is on the left axis, orange on the right. I'm using scipy's implementation of Runge-Kutta 45 in python to solve the ODE and have an array of $$A(t)$$ and $$B(t)$$ before I start (IE it's not a function, but derived from arbitrary sampled data).

I can get an accurate solution as long as I impose a small enough max_step parameter to catch the quick changes, but this takes a long time and I'm in a performance sensitive application. As an example, I need ~100k function evaluations for this function to have decent convergence (1E-4 accuracy). Of course the solver is wasting a huge amount of time/function calls in the regions where the coefficients aren't changing. I also have the arrays of coefficients before hand and in a sense already know where the trouble spots are and should be able to tell the solver about them before hand.

My question is: how do you efficiently solve an ODE with impulsive coefficients? Alternatively, is there a better way of solving the ODE on sampled data than using RK45? I simply interpolate the coefficients and pass it to the solver like a continuous function, but I'm not sure what others would recommend.

Edit: Another thought, is there a way to solve this ODE in quadrature kind of like how the solution to $$y^\prime(x)=f(x)y(x)$$ is just $$e^{\int f(x)dx}$$? I can't see an easy way even through diagonalizing $$M$$ because it depends on $$t$$, but who knows? Actually, does that equation hold for matrices w/ matrix exponentiation? It's not as simple as that, is it?

• Look into smoothing the problem like is shown in this blog post, and then use an ODE solver that's made for stiff equations. Also, setting the tstops or d_discontinuities at the jump points helps the solver hit the points exactly (though those options are not available in all software) Feb 5, 2020 at 5:27
• That's the key point: For accurate solutions, your time integrator has to put discrete time steps on all switch points. Feb 5, 2020 at 14:40
• Thanks! @WolfgangBangerth, is there a good way of translating my array of coefficients into a set of required time points? I saw something like this question which has an interesting solution to finding the minimum number of points required to represent the coefficients accurately with. Feb 5, 2020 at 22:40
• I think the first question to ask is where that set of coefficients come from. It's easier to think about why discontinuities are where they are than to try and reconstruct their locations at a later time. What is the physical mechanism that underlies these coefficients? Feb 6, 2020 at 0:04
• I'm modeling a physical system related to the transverse motion of relativistic charged particles in the presence of EM fields. The coefficients above are for a test system where the beam of particles encounters an accelerating field for the first ~0.2 ns, then a magnetic lens at ~0.7-1.2 ns. The spike is where the electric field ends and Maxwell's equations require a radial field to exist which gives a fast kick to the particles. For real systems I will use values from a numerical field solver which gives it to me at evenly spaced intervals along the problem domain, hence sampled data. Feb 6, 2020 at 2:51

I would recommend 3 things to pursue:

1. Use ode15s (a stiff solver). This will allow for better resolution near the jumps

2. Rescale the problem. Your coefficients are very large and if they are changing by several orders of magnitude, this is likely causing undue conditioning problems with your solver so rescaling then undoing the scaling after the ode solve is likely the best

3. If you are still getting bad errors at the jumps, try to specify events in ode15s options. You can specify events at the times of the jumps in the coefficients which will force the solver to solve near those points with high accuracy to resolve the discontinuity.

The phenomenon you're dealing with here is called a "stiff system". Check here for a full explanation: http://www.scholarpedia.org/article/Stiff_systems

Matlab’s ode45 is a non-stiff solver, thus ill-suited for your problem. You should try to use one stiff solver. Matlab has 4 stiff solvers but other tools like Scipy’s solve_ivp also possess non-stiff solvers as well: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

Julia also has solvers for stiff problems, for more details please refer to Chris Rackauckas' comment to your question above, since he is THE expert on the topic.

• Stiffness does create similar problems, but the OP's situation is actually much more harmless because the jumps don't arise from the system behaviour but are already fixed in time beforehand. So a stiff solver shouldn't be necessary; just refine the time step around the jump points. Feb 5, 2020 at 16:14