# Implicit integrator for ODE with quadratic right-hand side

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.

• If the quadratic term is not stiff, one possibility could be to perform a Lie or Strang splitting, retaining your backward Euler solution method for the linear term, and solving the evolution associated with the quadratic term via an explicit Runge-Kutta method. The downside is that you need to split your system (depends on how it is currently implemented), and that splitting methods are practically limited to second-order accuracy in time, but that may already be sufficient. Oct 7, 2021 at 9:32

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$$.