Implicit integrator for ODE with quadratic right-hand side

Question

I have an ODE for an unknown $x(t):[0,\infty)\to\mathbb R^n$ of the following form: $$ x_i'(t)=a_i^\top x(t) + x(t)^\top Q_i x(t), $$ for $i\in\{1,\ldots,n\}$. Here, the vectors $a_i\in\mathbb R^n$ form the columns of an $n\times n$ negative semidefinite matrix $A\in\mathbb R^{n\times n}$, and $Q_i\succeq0$ for all $i\in\{1,\ldots,n\}$. In words, the right-hand side of the ODE is a convex quadratic function of the state. We can write $x'=Ax+\mathcal Q(x)$, where each entry of $\mathcal Q$ is a convex homogeneous quadratic in $x(t)$.

Suppose $x'=Ax$ is well-approximated using backward (implicit) Euler integration to advance $t$, or even that we're willing to do exponential integration for this linear term. Is there a numerical integration method that might be effective in this case for timestepping the nonlinear ODE above?

As a simple example, integration based on variation of parameters could "bake in" solution of the linear part of the ODE. Or perhaps there is a variational version of implicit integration that would lead to a convex problem in our setting.

Answer 1

From your description, exponential integrators seem to be a good fit. They are based on variation of parameters and "bake in" an exact solution for the linear part of the ODE. There are many flavors of exponential methods that are determined by the discretization of the integral term in variation of parameters. Exponential Runge--Kutta methods tend to be the most popular. For example, exponential Euler for your problem is $$ x_{n+1} = x_n + h \varphi_1(h A) (A x_n + \mathcal{Q}(x_n)), $$ where $\varphi_1(z) = \frac{e^z -1}{z}$. Phi functions can be challenging to compute and are the core implementation issue. Fortunately, there are many resources available on this site and in the literature on this. If it is feasible to compute an eigendecomposition of $A$, I would strongly recommend exponential integrators.

I would not underestimate traditional methods like diagonally implicit Runge--Kutta or BDF methods. While they might not be as tailored to the problem, they are simpler, and I would still expect them to perform well. For the implicit stages, the nonlinear systems one must solve are quadratic, and it might be possible to use specialized solvers. I would also consider Rosenbrock-W methods such that you only need to solve linear systems with the SPD matrix $I - h \gamma A$.