A formal proof of the Lax equivalence theorem for finite difference schemes

Here $A$ is not assumed to be bounded, so that unbounded differential operators are included. The problem (1) is assumed to be well-posed, i.e., there exists a bounded, linear operator, $E\in B(Y,X)$, such that $EA=I$ in $D$, and that for $f\in Y$, equation (1) has a unique solution, $u=Ef$. Furthermore, the solution $u$ depends continuously on the data.

We assume that for each $h\in H$, problem (2) is well-posed and there exists a solution operator, $E_{h}=A_{h}^{-1}$, i.e. $u_{h}=E_{h}f_{h}$. The true solution $u$ and the approximate solution $u_{h}$ can be related with each other by defining a bounded, linear operator, $r_{h}:X\to X_{h}$ for each $h\in H$. Similarly, data $f\in Y$ can be related to data in a discrete space, $f_{h}\in Y_{h}$ by defining a restriction operator $s_{h}$. For each $h\in H$, $s_{h}:Y\to Y_{h}$ is also a bounded, linear operator. We assume that the operator norms can be uniformly bounded:

where the constants $C_{1},C_{2}$ are independent of $h$. The true solution $u=Ef$ is compared with the discrete solution $u_{h}=E_{h}s_{h}f$ corresponding to the discretized datum $f$. The family $(X_{h},Y_{h},A_{h},r_{h},s_{h})$ defines a method for the solution of (1) [33].

Intuitively, this means that in the limit of the discretization step, $h$, tending to zero, the numerical solution $E_{h}s_{h}f$ approaches the analytical solution $r_{h}Ef$. The analytical solution $r_{h}Ef$ is the restriction of the true (analytical) solution, $u=Ef$, onto the grid of size $N=1/h$, and $E_{h}s_{h}f$ is the discrete solution, $u_{h}=E_{h}f_{h}$ computed on the grid of size $N$.

Intuitively, this means that in the limit of the discretization step, $h$, tending to zero, the finite difference scheme $A_{h}u_{h}=f_{h}$ approaches the differential equation $Au=f$, i.e., we are discretizing the right differential equation.

Intuitively, stability of the numerical scheme means that a small numerical perturbation does not allow the solution to blow up. Uniform boundedness of the inverse $E_{h}=A_{h}^{-1}$ is a check on the conditioning of matrices (sensitivity to small perturbations), i.e., it ensures that the matrix $A_{h}$ is not ill-conditioned. Thus, if the numerical problem (2) were unstable, even though we were trying to solve the right differential equation, we would never converge to the true solution. Hence, both stability and consistency are sufficient for proving convergence of the numerical scheme.

The quantities within the norms (4) and (5) are, respectively, the global and local discretization errors.

In this Section we show how we formalized the proof of the Lax equivalence theorem [33] in the Coq proof assistant. All of the Coq formal proofs mentioned in this paper, containing the proofs of consistency, stability and convergence of finite difference schemes, and of the Lax equivalence theorem, are available at http://www-personal.umich.edu/~jeannin/papers/NFM21.zip.

The Coquelicot library [8, 7] defines mathematical structures required for implementing the proof. Since we use Coquelicot and standard reals library which are based on classical axiomatization of reals, our proofs are also non-constructive [8]. We define the Banach spaces (complete normed spaces, complete in the metric defined by the norm [25]) $(X,Y,X_{h},Y_{h})$ using a canonical structure, CompleteNormedModule, in Coq [16].

The definitions of the true problem (1) and the approximate problem (2) require that the mappings $A:X\to Y$ and $A_{h}:X_{h}\to Y_{h}$ be linear, and the solution operators $E:Y\to X$ and $E_{h}:Y_{h}\to X_{h}$ be linear and bounded. The linear mappings $A_{h}$ and $E_{h}$ are defined as functions of $h\in\mathbb{R}$. Boldo et al. [3] have defined linear mapping in the context of a ModuleSpace and bounded linear mapping in the context of a NormedModule in their formalization of the Lax Milgram Theorem in Coq [2, 15]. We extended these definitions in the context of CompleteNormedModule.

The definition of consistency (5) and convergence (4) hold in the limit of $h$ tending to zero. Thus, an important step in the proof is to express these limits in Coq. Formally, the notion of $f$ tending to $l$ at the limit point $x$ requires, for any $\epsilon>0$, to find a neighborhood $V$ of $x$ such that any point $u$ of $V$ satisfies $|f(u)-l|<\epsilon$ [8]. This notion has been formalized in Coquelicot [7] using the concept of filters. In topology, a filter is a set of sets, which is nonempty, upward closed, and closed under intersection [13]. It is commonly used to express the notion of convergence in topology. We have used a filter, locally x [27] to denote an open neighborhood of $x$, and predicate filterlim [27] to formalize the notion of convergence (in the context of limits) of $f$ towards $l$ at limit point $x$, i.e. $\lim_{x\to a}f(x)=l$. Therefore, the definition of consistency (5) is expressed as:

(is_lim (fun h:R => norm (minus (Ah h (rh h u)) (sh h (A u)))) 0 0

where the limits of functions is expressed using is_lim [8].

We next discuss the formalization of the statement of convergence of a finite difference scheme in Coq. We note that from Theorem 2.1, consistency and stability imply convergence. This notion is expressed in Coq as follows:

(is_lim (fun h:R => norm (minus (Ah h (rh h u)) (sh h (A u))))  0 0
    (*Consistency*) /\
(exists K:R , forall (h:R), operator_norm(Eh h)<=K ) (* Stability*) ->
is_lim(fun h:R=>norm (minus (rh h (E(f))) (Eh h (sh h (f))))) 0 0)
    (*Convergence*).

where the operator norm is defined as $||f||_{\phi}=sup_{u\neq 0_{E}\land\phi(u)}\frac{||f(u)||_{F}}{||u||_{E}}$ and has been formally defined in [3].

The basic idea is that we bound the global discretization error ($||r_{h}Ef-E_{h}s_{h}f||)$ above using the stability criterion, i.e. $||r_{h}Ef-E_{h}s_{h}f||\leq K||A_{h}r_{h}u-s_{h}Au||$, and then prove that as the local discretization error ($||A_{h}r_{h}u-s_{h}Au||)$ tends to zero in the limit of $h$ tending to zero, the upper bound on the global discretization error tends to zero (using the property of limits). Using the property of norm , i.e. $0\leq||r_{h}Ef-E_{h}s_{h}f||$, we arrive at the inequality

In Coq, we define the lower bound of the inequality as a constant function with value $0$ as: fun _ => 0. Since the limit of a constant function is the constant itself, i.e. $\lim_{h\to 0}0=0$, and $lim_{h\to 0}||A_{h}r_{h}u-s_{h}Au||=0$ (Consistency), using the sandwich theorem for limits, $\lim_{h\to 0}||r_{h}Ef-E_{h}s_{h}f||=0$. The sandwich theorem states that if we have functions obeying the inequality: $f(x)\leq g(x)\leq h(x)$ and $\lim_{x\to a}f(x)=L\quad\land\quad\lim_{x\to a}h(x)=L$ on some open neighborhood of $x=a$ , then $\lim_{x\to a}g(x)=L$. This proves the convergence of Definition 4 and completes the proof of the Lax equivalence theorem.

We want to prove that for a central difference approximation of the second derivative in the numerical scheme $\mathcal{N}_{h}$, the truncation error, $\tau=\mathcal{O}(\Delta x^{2})$. By invoking the definition of Big-O notation, the theorem statement can be stated as:

The equation (8) is stated formally in Coq as:

Theorem taylor_FD (x:R): Oab x ->exists gamma:R, gamma >0 /\ exists G:R,
G>0/\ forall dx:R, dx>0 -> Oab (x+dx) -> Oab (x-dx)->(dx< gamma ->
Rabs((D 0 (x+dx)- 2*(D 0 x) + D 0 (x-dx))*/(dx * dx)- D 2 x)<= G*(dx^2)).

where Oab x mean $a<x<b$ and D k x denotes $k^{th}$ derivative of $u$ with respect to x.
We start by introducing the following lemmas required to complete the proof.

Here, $F(x)$ is the Lagrange remainder in the expansion of $u(x+\Delta x)$ up to degree 3 and is defined as:

Thus, Lemma 1 states that the Lagrange remainder $F(x)=\frac{1}{4!}(\Delta x)^{4}\frac{d^{4}u(\xi)}{dx^{4}}$ is of order $(\Delta x)^{4}$ for all $\xi\in(x,x+\Delta x)$.

Here, $G(x)$ is the Lagrange remainder in the expansion of $u(x-\Delta x)$ up to degree 3 and is defined as:

Thus, Lemma 2 states that the Lagrange remainder $G(x)=\frac{1}{4!}(\Delta x)^{4}\frac{d^{4}u(\xi)}{dx^{4}}$ is of order $(\Delta x)^{4}$ for all $\xi\in(x-\Delta x,x)$.

Both the lemmas are a straightforward application of the Taylor–Lagrange theorem (Theorem 3.1), and are crucial to the formalization of the proof of the consistency of the finite difference scheme.

Next, we present an informal proof of the theorem followed by a discussion on the formal proof of the consistency theorem.

An important point to note is that the condition $|F(x)|+|G(x)|\leq M(\Delta x)^{4}+K(\Delta x)^{4}$ holds when $0<|\Delta x|<\gamma$, where $\gamma$ is as defined in (8). We therefore choose, $\gamma=min(\eta,\delta)$, where $\eta$ is such that, $|F(x)|\leq M(\Delta x)^{4}$ holds when $0<|\Delta x|<\eta$, and $\delta$ is such that, $|G(x)|\leq K(\Delta x)^{4}$ holds when $0<|\Delta x|<\delta$.

We followed the proof above and formalized it in the Coq proof assistant. To apply the Taylor–Lagrange theorem [28] to the consistency analysis of a central difference approximation, we broke down the theorem statement into two lemmas as discussed in the previous section. Therefore, in this section, we will discuss the proof of Lemma 1 and 2.

In this section, we relate the proof of consistency from Section 3.1 with the Lax equivalence Theorem 2.1. The numerical discretization of the differential equation can be expressed in the discrete domain as:

Comparing with the statement of consistency (5), we have

Taking the vector norm in the $L_{1}$ sense, $||.||_{1}$, equation (16) can be written as:

$\lim_{h\to 0}\frac{u_{o}}{h^{2}}=0$ and $\lim_{h\to 0}\frac{u_{N}}{h^{2}}=0$, trivially because of the boundary conditions we imposed, i.e. $u_{o}=0$ and $u_{N}=0$. The norm used in (16) are in the space $Y_{h}$, i.e., $||.||_{Y_{h}}$.
This reduces to proving:

But from the Taylor–Lagrange analysis discussed in section (3.1), we have

where $C$ is a constant, and $u_{i}=u(x_{i}),u_{i-1}=u(x_{i}-h),u_{i+1}=u(x_{i}+h)$. Substituting $\left.\frac{d^{2}u}{dx^{2}}\right|_{x_{i}}=1$, and using the inequality (19) and equation(18), we get

But, $\sum_{i=1}^{N-1}\lim_{h\to 0}|Ch^{2}|=0$. Hence, using the sandwich theorem, we prove that

In order to represent, $x_{i},\;i=0..N$, we define $x$ of type: nat $\to$ R. The boundary conditions are imposed as hypothesis statements:

    Hypothesis u_0 : (D 0 (x 0))= 0.
    Hypothesis u_N: (D 0 (x N)) =0.

The differential equation is defined as:

Hypothesis u_2x: forall i:nat, (D 2 (x i)) =1.

Equation (18) is formalized as a lemma statement:

Lemma lim_sum:is_lim (fun h:R =>
  sum_n_m (fun i:nat =>Rabs (( D 0 (x i -h) -2* (D 0 (x i))
    + D 0 (x i +h))*/(h^2) -1)) 1%nat (pred N)) 0 0.

This is where we integrate the proof of pointwise consistency of the FD scheme from section (3).

The main theorem statement which is an application of the statement of consistency required in the proof of Lax equivalence theorem from section (2) is as follows:

Theorem consistency_inst: forall (U:X) (f:Y) (h:R) (uh: Xh h)
 (rh: forall (h:R), X -> (Xh h)) (sh: forall (h:R), Y->(Yh h))
 (E: Y->X) (Eh:forall (h:R),(Yh h)->(Xh h)),
 is_lim (fun h:R => norm (minus (Ah h (rh h U)) (sh h (A U)))) 0 0.

We note here that the above-mentioned formalization is not unique to the second order scheme that we discussed. The approach we discuss can easily be generalized to verify consistency of any finite difference scheme. The crucial step in such a generalization is the appropriate instantiation of the $A_{h}$ matrix and the vectors $r_{h}u$ and $s_{h}Au$.

Analytical expressions for the eigenvalues and eigenvectors of $A_{h}(a,b,c)$ are given by:

$\forall m,j=1..N$. In Coq, we defined $\lambda_{m}$ and $s_{m}$ as follows:

Definition Eigen_vec (m N:nat) (a b c:R):= mk_matrix N 1%nat (fun i j =>
    sqrt ( 2 / INR (N+1))*(Rpower (a */c) (INR i +1 -1*/2))*
        sin(((INR i +1)*INR(m+1)*PI)*/INR (N+1))).

Definition Lambda (m N:nat) (a b c:R):= mk_matrix 1%nat 1%nat (fun i j =>
    b + 2* sqrt(a*c)* cos ( (INR (m+1) * PI)*/INR(N+1))).

Since naturals in Coq start with 0, we write INR (m+1) and INR i+1.

We then formally verify that the analytical expressions for the pair $(\lambda_{m},s_{m})$ indeed belong to the spectrum of $A_{h}$. From now on, we will refer to $A_{h}(a,b,a)$ as $A_{h}$ for the sake of brevity. In Coq, we state this formally as:

Lemma eigen_belongs (a b c:R): forall (m N:nat), (2 < N)%nat ->
    (0 <= m < N)%nat -> a=c /\ 0<c-> (LHS m N a b c) = (RHS m N a b c).

where, $LHS\overset{\Delta}{=}A_{h}s_{m}$ and $RHS\overset{\Delta}{=}s_{m}\lambda_{m}$. Here we used the definition of eigenvalue-eigenvector, i.e., $A_{h}s_{m}\overset{\Delta}{=}\lambda_{m}s_{m}$. Formalizing the proof of the lemma eigen_belongs was challenging due to the structure of the matrix $A_{h}$. $A_{h}$ is a tri-diagonal matrix with non-zero entries on the diagonal, sub-diagonal and super-diagonal. The other entries are zero and hence the matrix is sparse.

In Coq, we have to carefully destruct the matrix $A_{h}$ to separate the non-zero and zero sums in the LHS of equation (21). The idea is to do a case analysis on the row-index $i$, and has been illustrated in figure (1) in the Appendix 0.B.4. Details on the formal proof of the zero and non-zero cases are presented in Appendix 0.B.4.

Next, we discuss formalization of the boundedness of the matrix norm of $E_{h}=A_{h}^{-1}$. We have used an explicit formulation of $A_{h}^{-1}$ [17] in our formalization and we verify this formally using the definition: $A_{h}^{-1}A_{h}=I\;\land\;A_{h}A_{h}^{-1}=I$. Details on the proof can be referred to in the Appendix 0.B.1.

Here, we have used the definition of the spectral (2-norm): $||A||_{2}=\rho(A)$, where $\rho(A)$ is the spectral radius of $A$ and is defined as the maximum eigen-value of A, i.e. $\rho(A)=max_{m}|\lambda_{m}(A)|$. For the symmetric tri-diagonal matrix $A_{h}$, $A=E_{h}$ and $\lambda_{m}(E_{h})=1/\lambda_{m}(A_{h})$. Since $\lambda_{m}(A_{h})<0$, $max_{m}|\lambda_{m}(E_{h})|=1/|\lambda_{min}(A_{h})|$. Hence, we define the matrix norm in Coq as follows:

Definition matrix_norm (N:nat):= 1/ Rabs (Lambda_min N).

To show that the matrix norm is uniformly bounded, we need to show that $1/|\lambda_{min}(A_{h})|$ is uniformly bounded. This is where we instantiate the tri-diagonal matrix $A_{h}$ with the scheme $\mathcal{N}_{h}$. Thus, we prove the following lemma in Coq:

Lemma spectral: forall(N:nat),(2<N)%nat -> 1/Rabs(Lambda_min N) <= L^2/4.

where $L$ is the length of the domain, independent of $h$, and is constant throughout. Lambda_min is the minimum eigenvalue for the instantiated matrix, $A_{h}^{\prime}=A_{h}(\frac{1}{h^{2}},\frac{-2}{h^{2}},\frac{1}{h^{2}})$. We provide a paper proof of this bound in the Appendix 0.A.

To show that all the eigenvalues have the same bound, we prove that $\frac{1}{\lambda_{min}(A_{h}^{\prime})}$ is the maximum eigenvalue of $E_{h}^{\prime}$. The lemma statement is as follows:

Lemma eigen_relation: forall (i N:nat), (2<N)%nat ->(0<=i<N)%nat ->
    Rabs (lam i N) <= 1/ Rabs( Lambda_min N).

This completes the proof on the boundedness of the eigenvalues of $E_{h}^{\prime}$. The lemma, eigen_relation also shows that the spectral radius of $E_{h}^{\prime}$ is $\frac{1}{|\lambda_{min}(A_{h}^{\prime})|}$, and justifies the defintion of matrix_norm.

We note that the definition of the matrix norm of $A_{h}^{-1}$ is valid only if $A_{h}^{-1}$ is a normal matrix . We therefore verify that $A_{h}^{-1}$ is normal. The lemma statement is provided in the Appendix 0.B.3.

We also provide the proof that $A_{h}$ is diagonalizable in the Appendix 0.C. This helps us to formally establish that the eigen vectors are orthogonal and hence the eigen space is complete.

In this section, we integrate all of the previous lemmas to prove the main stability theorem (6).

Theorem stability: forall (u:X) (f:Y) (h:R) (uh: Xh h)
    (rh: forall (h:R), X -> (Xh h))(sh: forall (h:R), Y->(Yh h))
    (E: Y->X) (Eh:forall (h:R), (Yh h)->(Xh h)),
    exists K:R , forall (h:R), operator_norm(Eh h)<=K.

where the operator norm is instantiated with the matrix norm using the following hypothesis:

Hypothesis mat_op_norm: forall (u:X) (f:Y) (h:R) (uh: Xh h)
    (rh: forall (h:R), X -> (Xh h))(sh: forall (h:R), Y->(Yh h))
    (E: Y->X) (Eh:forall (h:R),(Yh h)->(Xh h)),
    operator_norm (Eh h) = matrix_norm m.

The total length of the Coq code and proofs is about 14,000 lines, of which about 1200 lines are specific to the scheme. The rest of the formalization can be reused for a generic symmetric tridiagonal matrix. It took us about 15 months for the entire formalization. Much of the effort was spent on destructing the matrices and developing required linear algebra tools to handle the matrix manipulation. Since we are treating stability from a spectral point of view, lack of spectral theory for numerical analysis for the Coquelicot definition of matrices has been challenging for us. For the proof of consistency, the primary challenge was the right placement of the quantifiers to bound the Lagrange remainder using the definition of big-$O$ notation. To instantiate $\Gamma=M+K$, we had to carefully destruct the lemmas into the main theorem. We believe that a generic library with an automated implementation of the big-O definitions would save considerable effort here. We also encountered issues in selecting appropriate instantiations for other existential parameters. In the proof of convergence, we had to carefully construct the application of properties of limit with filters of neighborhoods.