December 5, 2018 — Support my next blog post, buy me a coffee ☕.
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}. $$
2020 Exponential integrators for stiff PDEs
2018 Computer-assisted proofs for PDEs
2018 Spherical caps in cell polarization
2018 Solving nonlocal equations on the sphere
2018 Gibbs phenomenon and Cesàro mean
2017 Solving PDEs on the sphere
2017 When planets dance