# How to numerically solve noisy ODE

I need to numerically integrate a large number of ODE's of the following form $$\dot{X} = k_{1}\left[\rule{0pt}{4mm}U\left(t\right) - X\right] + k_{2}V\left(t\right)\left[\rule{0pt}{4mm}W\left(t\right) - X\right]$$

Here

• All variables and constants are positive real numbers
• $$X(t)$$ is the unknown function, that needs to be determined by integrating the ODE on the known finite interval $$t \in [t_0, t_1]$$.
• $$k_1$$ and $$k_2$$ are known positive constants
• $$U(t), V(t), W(t)$$ are known forcing functions, measured empirically on the interval $$t \in [t_0, t_1]$$ with a small uniform timestep $$\Delta t$$.

I don't need a very precise solution. A solution that has 5% error is perfectly tolerable, as long as the error does not explode towards the end of the integration interval.

I have used Scipy's SOLVE_IVP method to integrate the above equation, using different solvers, such as Runge-Kutta-45, and Radau. Since the above methods require knowledge of RHS of the ODE at arbitrary timesteps, I have used linear interpolation for the functions $$U(t), V(t), W(t)$$.

I have used numerical integration both on toy examples with smoothly changing forcing functions and known analytical solutions, as well as on real data. It works as expected for both. However, orders of magnitude more compute time is required to integrate real data.

I suspect that the main culprit is the noise in the measured forcing functions. On the one hand, it is undesirable that noise is not explicitly modeled by the solver: I suspect it could lead to systematic errors in the solution. On the other, noise in forcing functions likely results in adaptation mechanisms in the integrator (e.g. RK45) to pick very small step sizes, which result much slower computation.

I see 2 possible approaches to remedying this issue:

1. Low-pass filter the data, and only model the slowly-changing effects. It is given that interesting effects happen at slower timescale and noise at faster timescale, but I do not know the exact timescale cutoff, there is likely at least some overlap.
2. Include noise into the model, and use an integrator designed to handle noisy ODEs

Questions:

1. Is de-noising prior to integration at least somewhat feasible for this problem? Or do I have to explicitly model noise in my integration routine?
2. If yes, is there a rule of thumb to estimate the expected performance of a numerical integration scheme? E.g. how does the "smoothness" of empirical data affect the step size, chosen by the above integration routines. Would it help if I used a different interpolation routine (quadratic, spline, etc).
3. What methods exist that are designed for numerically integrating noisy ODEs?
• What do you mean by noisy? Could you show us an example? Commented Jun 20 at 11:12
• @Cesareo I cannot share closed-source data for legal reasons. I could try to synthesize data with similar power spectral density, if you wish. Or I could just find the PSD, and give you that. I'm just not sure how that would help you - could you please tell me? It is a complicated process that changes on multiple timescales. It has a clear slowly-changing component, and something that looks like "gaussian" noise on top of it. Judging by eye, the SNR is pretty good, noise is maybe 5-10% of the signal magnitude. Commented Jun 20 at 11:18
• $$X(t) = \exp \left(-\int _0^t(k_1+k_2 V(\xi ))d\xi \right) \int _0^t\exp \left(\int _0^{\eta }(k_1+k_2 V(\xi ))d\xi \right) (k_1 U(\eta )+k_2 V(\eta ) W(\eta ))d\eta +c_1 \exp \left(-\int _0^t(k_1+k_2 V(\xi ))d\xi \right)$$ Commented Jun 20 at 13:16
• No way, this is analytically solvable?? I'd never have imagined. Wasn't really the question in this post, but to me this is the best possible outcome. May I ask what technique you used? I guess it is a bit more complicated than basic integrating factors. Commented Jun 20 at 13:36

The ode

$$\dot X(t) +(k_1 + k_2V(t))X(t) = k_1 U(t)+k_2V(t)W(t)$$

is linear. It can be solved as the sum of the homogeneous solution $$X_h(t)$$ and a particular solution $$X_p(t)$$

$$\cases{ \dot X_h(t) +(k_1 + k_2V(t))X_h(t) = 0\\ \dot X_p(t) +(k_1 + k_2V(t))X_p(t) = k_1 U(t)+k_2V(t)W(t) }$$

The homogeneous solution is straightforward:

$$X_h(t) = c_0 \exp \left(-\int _0^t(k_1+k_2 V(\xi ))d\xi \right)$$

and now making $$X_p(t) = X_h(t)c_0(t)$$ we have

$$\dot c_0(t) = (k_1U(t)+k_2V(t)W(t))\exp \left(-\int _0^t(k_1+k_2 V(\xi ))d\xi \right)$$

and after integration

$$c_0(t) = \int _0^t(k_1 U(\eta )+k_2 V(\eta ) W(\eta )) \exp \left(\int _0^{\eta }(k_1+k_2 V(\xi ))d\xi \right)d\eta$$

hence $$X_p(t) = c_0(t)X_h(t)$$. Finally

$$X(t) = X_h(t) + X_p(t)$$

now when $$t = 0$$ we have $$X(0) = c_0$$

• I have checked it, all derivatives work out. There is a tiny discrepancy in the final answer. Since $X_p(t) = X_h(t)c_0(t)$ then $X(t) = X_h(t) + X_p(t) = X_h(t) (1 + c_0(t))$. Now, since $c_0(t)$ is an integral from $0$ to $t$, then $c_0(0) = 0$, and then $X(0)$ is just $c_0$. But this is really minor, just a constant definition Commented Jun 20 at 16:23
• In the homogeneous component, $c_0$ is a constant and is the ode constant to determine initial conditions. We used $c_0(t)$ as a way to determine $X_p$. We could instead call it as $\phi(t)$ to avoid confusion. Sorry. Corrected the value of $X(0)$. Commented Jun 20 at 16:57
• If I may, an important follow-up question. This is still a simplified version of my problem. In the current case, we are solving $\dot{X} + A(t)X = B(t)$, where everything is a scalar function. My ultimate problem has the same shape, but X is a 4d vector, and A,B are sparse 4x4 matrices. Can I hope to achieve a similar solution for the vector case? Just tell me yes or no, I'll figure it out if it is at all possible. I have seen that it is relatively simple if A and B are constant matrices, but if they are not I am not sure. Commented Jun 20 at 17:03
• The process to obtain the solution is the same. Commented Jun 20 at 17:20
• I have a bit of a problem. I have tried to numerically compute the value of the solution, and for reasonable bounds it seems to work fine. However, once I increase the integration time to a large enough value, the exponent in $c_0(t)$ is so big causes numerical overflow. The parameters are such that the solution should be finite. I suspect, there is some sort of singular behaviour, where the negative and positive exponent must cancel out. But I don't see at the moment what is the common constant multiplier that should cancel out. Commented Jun 21 at 19:52