Gabriel Rémond-Tiedrez

Dynamical systems is a mathematical field focused on studying the evolution of objects over time. For example, think of a particle of water in a river, a molecule being created in a chemical reaction happening in some medium, or the quantity of mature cells in blood. One usually starts with some rule for how to make a small step in the future, with the objective of understanding how these rules affect the very far future. Taking very small time steps is easy enough, we just need to follow our rules, the main difficulty is going beyond that. The good news is that Henri Poincaré existed — and of course many others before and after him — and Poincaré brought some very powerful and beautiful intuition linking dynamical systems with geometry. This allowed us to find answers to these questions in very elegant ways, and now serves as the basis for the study of dynamical systems. The idea is to really picture the time evolution as particles following some flow, and to consider the structures that this flow preserves. Studying the geometry of these structures and how they interact with each other then allows a complete understanding of the dynamical system. The following are 2 great online tools to play around with visualizing solutions to differential equations: a Phase Plane Plotter for Ordinary Differential Equations and VisualPDE for Partial Differential Equations.

One of the broad questions I am interested in is stability. A solution is called stable if solutions starting near it stay near it for all forward times, and is called asymptotically stable if solutions starting near it go towards it in forward time. The reason these solutions are relevant is because they are those observed through experiments, and proving the stability of a solution gives you important information about the behaviour of solutions at large times. One tool to show stability is a Lyapunov function, which one should think of as some decreasing energy, causing all solutions to go towards a stable state minimizing said energy. There are no general analytical tools to find such functions but thanks to Pablo Parrilo, we have a numerically tractable way of finding these in certain cases, using a field of convex optimization called Sum-Of-Squares (SOS) optimization. This is now a well established tool for obtaining numerical results, which my work aims to validate using Computer-Assisted methods, thus obtaining rigorous mathematical statements.

For the geometrically curious reader, Felix Klein introduced and solidified the idea that geometry is the study of invariant structures under some group of transformations acting on a space. In the case of continuous dynamical systems — differential equations — the flow of time really is the group \(\mathbb{R}\) acting on the space, where the invariant structures are the full orbits. This goes further into the more abstract study of dynamical systems, where one considers the action of non-compact groups on a state space. For example one could consider the action of \(\mathbb Z\) or \(\mathbb N\) on a space, which corresponds to a discrete dynamical system, where one applies a function over and over again to some space.

Some problems in dynamics can be reduced to showing polynomial inequalities. Presented below are 2 such classes of problems and an introduction to the Sum-Of-Squares methods used to deal with them, followed by an approach to validate these results, my main contribution.

Lyapunov functions: Some classical methods of showing stability are linear analysis or energy estimates. The former entails computing the linear part of the differential operator, be it finite or infinite dimensional, and controlling its spectrum in the right half-plane; the latter is a higher order method where one bounds energy of perturbations, usually in some \(L^2\) sense. However, these methods can be quite limited for interesting and complicated problems. Another option is looking for Lyapunov functions, which one can see as some modified energy. In particular, when looking at the ODE \(\dot x=f(x)\) in \( \mathbb R^n \) with an equilibrium at \(x=0\), a Lyapunov function is some function \( V:\mathbb R^n\to\mathbb R \) satisfying both inequalities \[ V(x)>0,\quad \frac{d}{dt}V(x(t))<0 \text{ for all } x\neq 0. \] While a method to find Lyapunov functions analytically does not exist for general differential equations, for polynomial vector fields, both conditions are reduced to a polynomial inequality.

Bounding maximal time-averages: In many chaotic systems, the behaviour of a single solution over some finite time is simultaneously very hard to obtain, and not particularly instructive. For the example of airflow over an airplane wing, one cares much more about the average lift over this wing than the lift at one specific point of a trajectory. As such, for some observable \(\varphi:\mathbb R^n\to\mathbb R\), one can consider the time-average starting at \(x_0\) as \[ \bar\varphi(x_0) = \limsup_{T\to\infty}\frac{1}{T} \int_0^T \varphi(x(t))\, dt, \text{ where } x(0)=x_0. \] Another advantage of focussing on averaged quantities is that transient behaviour gets eliminated. Moreover, one might also be interested in the extremal time-averages, that is \(\max_{x_0} \bar\varphi(x_0) \). Again taking an engineering perspective, this corresponds to a worst-case scenario, warranting careful investigation. Numerically computing these quantities directly requires integration over rather large times, a costly endeavour. However, with this wonderful article by Tobasco, Goluskin and Doering, the problem can again be tackled by reformulating this as a polynomial inequality, more precisely a convex polynomial optimization problem. This optimization procedure also gives localizing information about the maximizing trajectory, as can be seen in the video below.

SOS and SDP programming: The above results are dependent on our ability to show that \(p\ge 0\) for some polynomial \(p\). This is in general an NP-hard problem as the degree or the dimension of \(p\) grows, but one can replace positivity with the condition that \(p\) is a sum-of-squares (SOS). Indeed finding \(g_1,\dots,g_n\) such that \( p=\sum_{k=1}^n g_k^2 \) is a sufficient — but not necessary — condition for positivity. The real trick comes from observing that \(\sum_k g_k^2=z^T Q z\) for some vector of monomials \(z\) and positive semi-definite matrix \(Q\), and since the cone of positive semi-definite matrices is convex, this now becomes a convex feasibility problem: \(\exists Q\ s.t.\ Q\succeq0, A(Q)=b\). The linear equality constraints \(A(Q)=b\) correspond to matching coefficients in \(z^TQz\). To solve these semi-definite programs (SDP), one uses interior-point method algorithms, which can give good numerical results in polynomial time. However, the outputs usually violate the equality constraints by some numerical tolerance, meaning a layer of post-processing is required to rigorously establish positivity.

Validating results: Upon inputting a polynomial \(p\) and a choice of monomials \(z\), SDP solvers will return a positive definite matrix \(Q\succ 0\) such that \(p=z^TQz+r\), where \(r\) is small. This \(r\) is some remainder polynomial with size corresponding to the numerical tolerance allowed by the solver, and this is the crux of the problem. Indeed \(Q\succ 0\) only gives \(p-r\ge0\), but \(r\) could very well be negative somewhere in the domain. Supposing we can write \(r=z^TRz\) for some symmetric error matrix \(R\), the problem is reduced to showing that the sum \(Q+R\) is positive semi-definite. In principle, one expects this to hold only if \( \lambda_{\min}(Q)> \|R\|_2 \), so that \(Q\) is positive enough that the small error \(R\) is not enough to push the eigenvalues below zero. We therefore have this interplay between the smallest eigenvalue of \(Q\) and the size of the numerical error. One difficulty is that solutions usually lie very close to the boundary of positivity, moreover decreasing the numerical tolerance usually also decreases the smallest eigenvalue, making it very difficult to obtain numerical results apt for validation. In fact, one can find examples of relatively low-dimensional polynomials which numerically appear to be sums-of-squares, but for which validation appears out of reach. This failure of validation emphasises the dependence on high quality numerical results, the main obstacle here. The good news is that when looking at examples from dynamics, one can modify parameters in order to get solutions further away from the boundary, thus allowing validation.