Solving numerically a linear ODE

Question

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: \begin{equation} \label{eq:method} y'[n] \simeq (\alpha+1) f_s \ y[n] - (\alpha+1) f_s \ y[n-1] - \alpha \ y'[n-1], \end{equation} 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?

Thank you in advance.

Answer 1

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/