December 5, 2018
A couple of weeks ago, I spent a few days at McGill University with Jean-Philippe Lessard. He introduced me to the world of computer-assisted proofs, which I briefly describe in this post.
Suppose we want to find steady-state solutions of nonlinear time-dependent PDEs of the form $$ u_t = F(u), $$ for some smooth operator \(F:U\mapsto U'\) between two Banach spaces \(U\) and \(U'\). This amounts to solve the nonlinear equation \(F(u) = 0\). We typically use a spectral method and interpret \(u\) as an infinite vector of spectral coefficients and \(F\) as an operator acting on those coefficients.
For finite \(N\)-dimensional problems \(F_N(u_N)=0\) one would typically use Newton's method where one iterates the map $$ T_N(u_N) = u_N - DF_N^{-1}(u_N)F_N(u_N). $$ When the Jacobian \(DF_N^{-1}(u_N)\) is invertible, \(T_N\) is a contraction and it allows us to prove that the numerical solution \(\bar{u}_N\) is close to the true solution \(u_N\). (Typically, \(u_N\) would be the vector of coefficients in a spectral expansion.)
In an infinite-dimensional setting, Newton's method reads $$ T(u) = u - DF^{-1}(u)F(u). $$ However, working with the inverse \(DF^{-1}(u)\) of the Fréchet derivative might be hard. Instead, one typically considers an approximation \(A\) (independent of \(u\)) of \(DF^{-1}(u)\) and looks at the fixed point problem $$ T(u) = u - AF(u). $$ The main challenge is to choose \(A\) such that \(T\) is a contraction on some neighborhood of the unknown true solution.
The idea behind computer-assisted proofs goes like this. First, compute a numerical solution \(\bar{u}_N\) of the discretized problem \(F_N(u_N)=0\) in \(N\) dimensions. Second, reinterpret \(\bar{u}_N\) as an element \(\bar{u}\) of the infinite-dimensional space by padding the vector of coefficients with an infinite number of zero. We expect the true solution \(u\) to be close to \(\bar{u}\) so we will choose \(A\) to be an approximation of \(DF^{-1}(\bar{u})\). The following theorem justifies this approach.
Theorem (Radii polynomial approach in Banach spaces). Let \(U\) and \(U'\) be Banach spaces. Denote the norm on \(U\) by \(\Vert\cdot\Vert_U\). Consider bounded linear operators \(A^{\dagger}\in L(U,U')\) and \(A\in L(U',U)\). Assume that \(F:U\mapsto U'\) is \(C^1\), \(A\) is injective and \(AF:U\mapsto U\). Consider an approximate solution \(\bar{u}\) of \(F(u)=0\) (usually obtained using Newton's method on a finite-dimensional discretization). Let \(Y_0\), \(Z_0\), \(Z_1\) and \(Z_2\) be positive constants satisfying $$ \begin{align*} \Vert AF(\bar{u})\Vert_U & \leq Y_0, \\[1.5pt] \Vert I - AA^\dagger\Vert_{B(U)} & \leq Z_0, \\[1.5pt] \Vert A[DF(\bar{u})-A^\dagger]\Vert_{B(U)} & \leq Z_1, \\[1.5pt] \Vert A[DF(u)-DF(\bar{u})]\Vert_{B(U)}& \leq Z_2r, \end{align*} $$ for all \(\Vert u - \bar{u}\Vert_U\leq r\), where \(\Vert\cdot\Vert_{B(U)}\) is the induced operator norm for bounded linear operators from \(U\) to itself. Define the radii polynomial $$ p(r) = Z_2 r^2 - (1-Z_1-Z_0)r + Y_0. $$ If there exists \(r_0>0\) such that \(p(r_0)<0\), then there exists \(u\in B_{r_0}(\bar{u})\) satisfying \(F(u)=0\).
One of the main challenges in this framework is to construct such an operator \(A\). (Note that \(A^\dagger\) is an approximation of the Fréchet derivative \(DF(\bar{u})\). In practice we often start by selecting a good \(A^\dagger\).)
Consider the Allen–Cahn equation with periodic boundary conditions, $$ u_t = F(u) = \epsilon^2 u_{xx} + u^3 - u. $$ In Fourier space, the linear part \(Lu=\epsilon^2 u_{xx}-u\) is diagonal with entries \(-\epsilon^2k^2 - 1\), \(k\in\mathbb{Z}\).
Let us denote by \(\bar{u}_{N}\) the numerical solution obtained using wavenumbers \(k\) between \(-N\) and \(N\) (\(2N+1\) grid points) and by \(DF_N(\bar{u}_N)\) the \((2N+1) \times (2N+1)\) Jacobian matrix at \(\bar{u}_N\), which we can compute. In this case, a good choice for \(A^\dagger\) is $$ A^\dagger = \begin{pmatrix} \scriptsize{DF_N(\bar{u}_N)} & 0 & 0 & \ldots \\ 0 & \scriptsize{-\epsilon^2(N+1)^2 -1} & 0 & \ldots \\ 0 & 0 & \scriptsize{-\epsilon^2(N+2)^2 -1} & \ldots \\ \vdots & \vdots & \vdots & \ddots \end{pmatrix}. $$ The idea is that for large enough \(N\) the Fréchet derivative of \(L\) will dominate the Fréchet derivative of \(F\). (Note that \(DL(u)=L\).) Therefore we can choose $$ A = \begin{pmatrix} \scriptsize{DF_N^{-1}(\bar{u}_N)} & 0 & 0 & \ldots \\ 0 & \frac{1}{-\epsilon^2(N+1)^2-1} & 0 & \ldots \\ 0 & 0 & \frac{1}{-\epsilon^2(N+2)^2-1} & \ldots \\ \vdots & \vdots & \vdots & \ddots \end{pmatrix}. $$