# Solving numerically a linear ODE

I start by saying that I do not have a strong background in numerical analysis, so I may miss some basic things or make trivial mistakes.

Motivated by some problems in digital signal processing, I would like to solve numerically a linear ODE of order $$m$$ of the form: $$\begin{equation*} y(t)+ a_1(t) y'(t) + \dots + a_m(t) y^{(m)}(t) = b_0(t) x(t) + b_1(t) x'(t) + \dots + b_m(t) x^{(m)}(t), \end{equation*}$$ where $$x$$ is known (input) and $$y$$ is unknown (output). Let $$f_s > 0$$ be the sampling rate and $$y[n] = y(n/f_s)$$ the approximate solution we want. We can also put: $$\begin{equation*} f(t) := b_0(t) x(t) + \dots + b_m(t) x^{(m)}(t). \end{equation*}$$

This is what I have done so far. First of all, I reduced the equation to a linear first-order differential system in the usual way: $$\begin{equation*} \begin{cases} y_0(t) + a_1(t) y_1(t) + \dots + a_m(t) y_m(t) = f(t) \\ y'_0 (t) = y_1(t) \\ \dots \\ y'_{m-1}(t) = y_m(t). \\ \end{cases} \end{equation*}$$ Then I used the following linear 1-step method: $$$$\label{eq:method} y'[n] \simeq (\alpha+1) f_s \ y[n] - (\alpha+1) f_s \ y[n-1] - \alpha \ y'[n-1],$$$$ where $$\alpha \in [0,1]$$. So I got the following linear system $$\begin{cases} y_0[n] + a_1[n] y_1[n] + \dots + a_m[n] y_m[n] = f[n] \\ y_{i+1}[n] = (\alpha+1) f_s \ y_i[n] - (\alpha+1) f_s \ y_i[n-1] - \alpha \ y_{i+1}[n-1], \ i = 0, \dots, m-1 \\ \end{cases}$$ and the solution is recursively given by: $$\begin{cases} y_0[n] = \dfrac{f[n] + \sum_{i=1}^m \{ a_i[n] \sum_{j=1}^i (\alpha+1)^{i-j} f_s^{i-j}((\alpha+1) f_s \ y_i[n-1] + \alpha \ y_{i+1}[n-1])\}}{1 + \sum_{i=1}^m a_i[n] (\alpha+1)^i f_s^i} \\ y_{i+1}[n] = (\alpha+1) f_s \ y_i[n] - (\alpha+1) f_s \ y_i[n-1] - \alpha \ y_{i+1}[n-1], \ i = 0, \dots, m-1 \end{cases}$$

Here are my questions.

1) Did I do anything wrong?

2) Does the reduction from one $$m$$-th order equation to a system of $$1$$-st order equations affect the quality of the solution? Is there any alternative? For example, how do you compare it to multistep methods?

3) The most general linear 1-step method is given by: $$\begin{equation*} y'[n] \simeq k_1 y[n] + k_2 y[n-1] + k_3 y'[n-1], \end{equation*}$$ where $$k_1, k_2, k_3 \in \mathbb R$$. The previous expression in terms of $$\alpha$$ and $$f_s$$ is obtained asking the method to be $$1$$-st order at least, while the additional condition $$\alpha \in [0,1]$$ is used to avoid stability issues. Notice that for $$\alpha = 0$$ you get backward Euler, for $$\alpha = 1$$ the trapezoidal rule and for $$\alpha \to \infty$$ forward Euler. I found this class of methods in previous papers (e.g. Germain, François G., and Kurt J. Werner. "Design principles for lumped model discretisation using Möbius transforms." Proc. DAFx-15, 2015), but I can not find more information. Do you know any reference on that? Can you compare it, in terms of efficiency and accuracy, with more usual (classes of) methods?

4) In general, could you provide me references for numerical methods for ODEs, especially linear?

• There is so much literature on ODE solvers, including for free on Wikipedia. Have you looked through some of it? – Wolfgang Bangerth Apr 20 '20 at 20:41
• @WolfgangBangerth I read something, but the literature is so vast that I got a bit lost. I also didn't find the method I was talking about, that's why I thought that probably somebody more experienced could help. – avril_14th Apr 21 '20 at 6:30

It seems like your goal is to get an accurate numerical solution for you differential equation, which likely does not require you to code your own ODE solver. In that case, it is likely more efficient for you to reframe your problem for use in an ODE solver in your programming environment of choice.

A standard first step is to rewrite your equation as a first-order vector differential equation. Let $$v = (y,y',y'',\dots,y^{(m-1)})^T$$. We then have the system \begin{aligned} v_1' &= v_2 \\ v_2' &= v_3 \\ & \ \ \vdots \\ v_m' &= \frac{1}{a_m}\left(-(v_1+a_1v_2+a_2v_3+\dots a_{m-1}v_m)+g(t)\right) , \end{aligned} where $$g(t):=(b_0x+b_1x'+\dots+b_mx^{(m)})$$. I have suppressed dependence on $$t$$ for brevity.

There are a few issues/things to think about before plugging into an ODE solver.

1. You need an initial condition for each component of $$v$$, i.e., $$y$$ and its first $$m-1$$ derivatives. This is required not only for numerical algorithms but also for uniqueness of the true underlying solution.

2. What if $$a_m(t)$$ is $$0$$/is very very small? In this case, division by $$a_m$$ may not be valid so you will need to keep it on the lefthand side multiplying $$v_m'$$. This can still be solved, but you will need to use a solver capable of using time-dependent mass matrices. In this case, the mass matrix will be a diagonal matrix with $$m-1$$ ones and $$a_m(t)$$ for the last element.

3. How do I compute $$a$$'s and $$b$$'s at times for which I do not have data? Since you are coming from signal processing, I am assuming you only have access to these quantities at a discrete set of points. You can avoid this by using an ODE solver that uses a fixed time step and match this up exactly with your sampling points, but for anything with variable time stepping you will have to use some sort of interpolation. Depending on how finely your data was sampled, you may be able to get away with a simple linear interpolation. If more accuracy is needed, you can try cubic splines or something similar. If you do not have much background in numerical analysis, I would recommend using a pre-existing package for this as numerically stable interpolation is often tricky business

4. What about $$x$$? This will likely be the hardest part of this problem. Assuming it is determined from sampled data, you will have the same problems as you did with $$a$$ and $$b$$, i.e., interpolation may be necessary. However, since you need $$m$$-th derivatives of $$x$$, you need very smooth interpolants to ensure that these derivatives exist. For instance, you may need to use splines with order $$m+1$$ polynomials. For both this and interpolation of $$a$$ and $$b$$, you will likely have a problem with exrapolation outside of your sampled data, but there is likely nothing you can do about that.

If your data is not sampled, then ignore points 3 and 4. To answer your more general question, there are tons of numerical methods for solving ODEs and the correct choice of algorithm depends on particular pathologies of your particular problem. First order methods are rarely used except in very specific situations. You reference a paper who is using these methods, so there may be some domain-specific reason to do so, but the problem as you have presented it here does not seem to me like it warrants a first order scheme and higher order schemes should work just fine once you work out the details above.

This is written for the Julia package DifferentialEquations.jl, but many of the methods mentioned are implemented in ODE packages in other languages. Try reading this to get a sense of what solver may work best for your problem: https://diffeq.sciml.ai/latest/solvers/ode_solve/

• For the linear specific solvers, it may be useful to link over to the tooling library ExponentialUtilities.jl which is where all of the Krylov exponential methods are implemented. Our LinearExponential is a thin wrapper over things like KIOPS which is likely what they want. – Chris Rackauckas Apr 21 '20 at 17:42
• Their problem has time-varying coefficients so wouldn't something like a Magnus integrator be more appropriate than than exponential-based methods? Obviously you know the capabilities of the integrators in Julia better than anybody but I wasn't aware that LinearExponential could handle time-varying operators. – whpowell96 Apr 21 '20 at 20:37
• Oh, I missed that it was time-varying. Yes, with time-varying operators a Magnus method is more appropriate. You'll have order loss if you use LinearExponential on that. – Chris Rackauckas Apr 21 '20 at 21:09
• Thanks a lot. I was already aware of 1. and 2. Regarding 3. and 4., the common solution is using a fixed-step method where the step is exactly the sample rate, so that no interpolation is needed. In order to compute higher derivatives of x, I just wanted to use repeatedly the approximation formula above. Does it work or is there too much propagation error? – avril_14th Apr 22 '20 at 6:35
• I do not think that such a method would work. Generally, to approximate the $m$-th derivative, you need information from at least $m+1$ points. More importantly is that you want to ensure that the approximation of the derivatives of $x$ is more accurate than the ODE solver to ensure the error propogation is not out of control. This is especially important at the beginning of the solve, as you will need values of $x$ and all its derivatives to at least 2nd order at points where you have very few sampled points in the past – whpowell96 Apr 22 '20 at 14:46