System of ordinary differential equations - time complexity of initial value problem

Question

I am interested in knowing what the time complexity is (in Big-$\mathcal O$ notation) for solving system of $N$ differential equations?

I am using ode15s in MATLAB for this but have no clue about time complexity.

Answer 1

I'm going to assume that you're specifically talking about Initial Value Problems (IVPs) for Ordinary Differential Equations (ODEs), since that's what ode15s() is used for.

First, a disclaimer: Good ODE solvers have stepsize control procedures that automatically take smaller or larger steps to achieve a desired accuracy, which means that the actual amount of time taken to solve an IVP can vary a lot depending on the details of the problem. Also, the big-O notation hides constants which are often important when comparing performance of different methods.

That said, we can easily comment on the time complexity of each step taken by a solver. Supposing we have an IVP $$ \dot{x} = f(x,t), \qquad x(0) = x_0 $$ with $ x \in \mathbb{R}^N $, then for explicit solvers like Matlab's ode45(), the complexity per step is the same as the complexity of evaluating $f(x,t)$. This is because explicit solvers evaluate $f(x,t)$ a fixed number of times per step. Typically this means a complexity of $O(N^2)$, since each component of $f(x,t)$ can depend on all components of $x$, but is less in some cases. For example, if $f(x,t)$ is 'sparse' in the sense that each component only depends on at most $m$ components of $x$, then the time complexity per step drops to $O(mN)$.

On the other hand, if you use an implicit solver like ode15s() then a linear system of equations is solved at each step. This involves an $N \times N$ matrix, and the complexity of solving them is $O(N^3)$. This is asymptotically greater than the time required to evaluate $f(x,t)$, so is also the time complexity of the step. Once again though, sparsity can sometimes reduce this for a specific set of ODEs. If each component of $f(x,t)$ only depends on at most $m$ components of $x$, we can often construct the linear system such that it involves a banded matrix of bandwidth $m$, and the time complexity per step reduces to $O(m^2N)$.

Now here's where it gets complicated. If we were to fix a stepsize $\delta t$ and take $M$ steps of that fixed size to solve a particular IVP, then we could get the time complexity of the whole algorithm by multiplying the above step complexities by $M$, for example $O(MN^3)$ in the implicit case. But usually solvers take steps $\delta t$ of different sizes to achieve a desired accuracy. This makes it impossible to say anything about the overall time complexity without also considering details of the specific IVP. Despite the fact that implicit methods like ode15s() have a worse time complexity per step than explicit methods like ode45(), for some problems (called stiff problems) we find that implicit methods outperform explicit methods as they can take larger steps while maintaining accuracy. Attempts to analyse this usually consider the 'stability' of the numerical method used on a particular ODE for a particular stepsize, as well as the 'order of accuracy' of the method, and there's a very large literature on the subject.

In practice, the usual procedure is to try an explicit solver like ode45() (which is the Dormand-Prince method) and then if the problem is taking a long time to solve, try an implicit method like ode15s(). There's a huge number of different methods and software for solving IVPs, so usually the only way to find out the best method for a particular problem is to benchmark.