Formal verification of iterative convergence for differential equations]Formal verification of iterative convergence for numerical solutions of differential equations

[1]Mohit<EMAIL_ADDRESS>
1]Joshua<EMAIL_ADDRESS>
1]Jean-Baptiste<EMAIL_ADDRESS>
[1]Department of Aerospace engineering, University of Michigan,  Ann Arbor, 48109, MI, USA

Physical systems are usually modeled by differential equations, but solving these differential equations analytically is often intractable. Instead, the differential equations can be solved numerically by discretization in a finite computational domain. The discretized equation is reduced to a large linear system, whose solution is typically found using an iterative solver. We start with an initial guess, $x_0$, and iterate the algorithm to obtain a sequence of solution vectors, $x_m$. The iterative algorithm is said to converge to the exact solution $x$ of the linear system if and only if $x_m$ converges to $x$.

It is important that we formally guarantee the convergence of iterative algorithm, since these algorithms are used in simulations for design of safety critical systems such as airplanes, cars, or nuclear power plants.

In this paper, we first formalize the necessary and sufficient conditions for iterative convergence in the Coq proof assistant. We then extend this result to two classical iterative methods: Gauss-Seidel iteration and Jacobi iteration. We formalize conditions for the convergence of the Gauss–Seidel classical iterative method, based on positive definiteness of the iterative matrix. We then use these conditions and the main proof of iterative convergence to prove convergence of the Gauss-Seidel method. We also apply the main theorem of iterative convergence on an example of the Jacobi classical iterative method to prove its convergence.

We leverage recent developments of the Coq linear algebra, 's real analysis library and the  library for our formalization.

§ INTRODUCTION

Differential equations are used throughout science and engineering to model and predict the behavior and performance of aircraft, cars, trains, robots, chemical powerplants, the stock exchange, and many more.
For example, the aerodynamics of an airplane can be represented by the Navier–Stokes equations [1]; the Riccati differential equation [2] is encountered in the problems of optimal control; and the Black–Scholes equation [3] is a partial differential equation that is used to model valuation of stock options [4] in mathematical finance.
Differential equations are thus pervasive in almost every aspect of science and engineering, and being able to solve those differential equations precisely and accurately, but also while trusting that the solutions are accurate, is of utmost importance.

However, in many cases, solving the differential equations analytically is intractable.
Therefore, the differential equations are solved numerically in a finite computational domain, using methods such as finite differences or finite elements.
In many cases, the desired precision is so high that it requires a large finite domain, and the exact numerical computation – typically involving the inversion of a large matrix – cannot be performed exactly. In those cases, a common approach is to use a family of methods called iterative methods [5]. Instead of inverting a large matrix, those methods build a series of approximations of the solution that eventually converge to the direct solution.
Many general purpose ordinary differential equation (ODE) solvers use some kind of iterative method to solve the linear system. For instance, ODEPACK [6], which is a collection of FORTRAN solvers for initial value problem for ODEs, uses iterative (preconditioned Krylov) methods instead of direct methods for solving linear systems.

Since most codes in scientific computing are still written and maintained in FORTRAN, these solvers are being widely used. Another widely used suite of ODE solvers is SUNDIALS [7].
SUNDIALS has support for a variety of direct and Krylov iterative methods for solving the system of linear equations.

SUNDIALS solvers are used by the mixed finite element (MFEM) package for solving nonlinear algebraic systems and by NASA for spacecraft trajectory simulation [7].

Because those iterative methods are widely used, it is important to obtain formal guarantees for the convergence of iterative solutions to the “true” solutions of differential equations.
In this paper we use the Coq theorem prover to formalize the convergence guarantees for iterative methods.

A number of works have recently emerged in the area of formalization of numerical analysis. This has been facilitated by advancements in automatic and interactive theorem proving [8, 9, 10, 11]. Some notable works in the formalization of numerical analysis are the formalization of Kantorovich theorem for proving convergence properties of the Newton methods in Coq by Ioana Pasca [12], the formalization of the matrix canonical forms in Coq by Cano et al. [13], and the formalization of the Perron-Frobenius theorem in Isabelle/HOL [14] for determining the growth rate of $A^n$ for small matrices $A$. Brehard [15], worked on rigorous numerics that aims at providing certified representations for solutions of various problems, notably in functional analysis. Boldo and co-workers [16, 17, 18] have made important contributions to formal verification of finite difference schemes. They proved consistency, stability and convergence of a second-order centered scheme for the wave equation. Tekriwal et al. [19] formalized the Lax equivalence theorem for finite difference schemes. The Lax equivalence theorem provides formal guarantees for convergence of a finite difference scheme if it is stable and consistent. Besides Coq, numerical analysis of ordinary differential equations has also been done in Isabelle/ HOL [20]. Immler et al. [21, 22, 23] formalized flows, Poincaré map of dynamical systems, and verified rigorous bounds on numerical algorithms in Isabelle/HOL. In [24], Immler formalized a functional algorithm that computes enclosures of solutions of ODEs in Isabelle/HOL.
However, the problem of iterative convergence has not been dealt with formally before. Over the recent years, the  library [25] in Coq has made significant contributions in formalizing the algebraic structures and theorems and lemmas from linear algebra. The time is right to exploit and leverage those developments to address problems in verified scientific computing. We believe that our formalization on iterative convergence is an important step towards making fully verified scientific computing possible.

Overall this paper makes the following contributions:

* We provide a formalization of the necessary and sufficient conditions for iterative convergence in the Coq proof assistant;

* We formalize conditions for convergence of the Gauss–Seidel classical iterative method;

* We then apply the main theorem on iterative convergence to an example of the Jacobi iteration, another classical iterative method, to prove its convergence;

* During our formalization, we develop libraries for dealing with complex matrices and vectors on top of the  complex formalization;

* We also formalize the properties of the 2-norm of a matrix and its spectral properties.

This paper is structured as follows. In Section <ref>,
we discuss an overview of the existing libraries in Coq that we base our formalization on. In particular, we discuss those formalizations that are relevant to follow our work. We then define the modulus of a complex number, matrix and vector norms, which we use extensively in our formalization of iterative convergence. In Section <ref>, we formalize the necessary and sufficient conditions for the iterative convergence in the Coq proof assistant. We formally prove that an iterative method converges if and only if all the eigenvalues of an iterative matrix have magnitude less than 1. We then state an example problem in Section <ref> that we later use to demonstrate the application of the iterative convergence theorem to a couple of classical iterative methods. In Section <ref>, we formalize the sufficient conditions for convergence of the Gauss–Seidel iterative method. The sufficient condition states that all the eigenvalues of the iterative matrix of a real and symmetric matrix $A$ have a magnitude less than 1 if $A$ is positive definite. We then apply these conditions and the theorem from Section <ref> to formally prove that the sequence of iterative solutions $x_m$ obtained by Gauss–Seidel method converges to the solution $x$ obtained using direct matrix inversion methods. In Section <ref>, we apply the theorem on iterative convergence from Section <ref> to prove convergence of an example of Jacobi iterative method.

Both Gauss–Seidel and Jacobi iteration are instances of iterative methods. In Section <ref>, we discuss existing gaps in the current state-of-the art and provide a detailed account on how we addressed them in our formalization. In Section <ref>, we conclude by summarizing key takeaways and discussing future work.
Our full formalization is available at: <https://github.com/mohittkr/iterative_convergence.git>.

§ PRELIMINARIES

In this section, we discuss some preliminaries required to follow our work on formalizing iterative convergence.

§.§ Overview of the libraries in Coq

In our formalization, we use use the Standard reals library, the real analysis library from  [9], and the linear algebra package from  [25].
Before we discuss their use in our formalization, it is worth spending some time discussing important features from these libraries.

[9] is an extension of the standard real analysis library and formalizes important results in real analysis like limits, derivatives, integrals etc. In particular, we used the limits formalization from  [9] extensively. Limits of a sequence are defined in  [9] using the  predicate.
This predicate is implemented using filters and hence the definition of  is expressed in Coq as
Definition is_lim_seq (u : nat → R) (l : Rbar) :=
filterlim u eventually (Rbar_locally l).

There exists an equivalent $\epsilon - \delta$ definition of limits in  [9] which states that sequence $(u_n)$ has limit $l$ if and only if

\begin{equation*}
\forall \epsilon>0, \exists N \in \mathbb{N}, \forall n \in \mathbb{N}, n \geq N \Rightarrow |u_n - l| < \epsilon
\end{equation*}

Formally, this is stated in  [9] as
Definition is_lim_seq' (u : nat → R) (l : Rbar) :=
match l with
| Finite l => forall eps : posreal,
eventually (fun n => Rabs (u n - l) < eps)
| p_infty => forall M : R, eventually (fun n => M < u n)
| m_infty => forall M : R, eventually (fun n => u n < M)

where  is a topology on natural numbers and states that given a property $P$, there exists a natural number $N$ such that a property $P\; n$, holds for all $n \geq N$.
In Coq,
Definition eventually (P : nat → Prop) :=
exists  N : nat, forall n, (N <= n)

In our formalization, we often had to switch between the filter definition of limits and the $\epsilon-\delta$ definition of limits. This was possible due to the following equivalence lemma in  [9]
Lemma is_lim_seq_spec :
forall u l, is_lim_seq' u l <-> is_lim_seq u l.

We also work extensively with matrices in our formalization. Thanks to the  [25] formalization of linear algebra, we have all the tools to work with matrices. In  [25], matrices are defined as a function from finite sets to an appropriate ring type.
Variant matrix : predArgType := Matrix of ffun 'I_m × 'I_n → R.

Here  and  are ordinal types and represent a finite set of naturals from $\{0,.., m-1\}$ and $\{0,.., n-1\}$, respectively. The following notation lets us define a matrix as a function of  and .
Notation "_ ( i , j ) E" := (matrix_of_fun (fun i j => E)).

For instance, a $2 \times 2$ real valued matrix, $A = [1, 2; 3, 4]$ can be defined as
Definition A := _(i< 2, j < 2)
(if (i = = 0
(if (j = = 0
(if (j = = 0

Another formalization in  [25] that was really helpful to us was the  library for definite iterated operations. The notation
Notation "[ op / idx ]_ ( i <- r  P ) F" :=
(bigop idx r (fun i => BigBody i op P

allows us to define iterated sums and products by instantiating the  operator and the appropriate identity . Here,  is a function of  chosen from a finite sequence  when the predicate  holds true. The matrix operations like matrix-vector multiplication, dot products, traces etc. are defined in terms of these big operations.

The  library in  [25] allows us to define a finite sequence. In our formalization, we use sequences to get the eigenvalues and its multiplicity from the roots of the characteristic polynomial of a matrix A. Therefore, it is worth going through some relevant operations from the sequence library.
The following notation defines a map for each element  in the sequence .
[seq E | x <- s] := map (fun x => E) s.

To extract an $n^{th}$ element in the sequence, we use the following notation from
nth x0 s i

The keyword  allows us to translate from the set notation to sequences.

Another notable formalization that we used in our formalization was the canonical forms library [13]. The canonical forms library defines Jordan canonical forms which we use extensively in our proof of iterative convergence. In particular, we will be using the definition of a block diagonal matrix which is defined in Coq as:
Definition diag_block_mx s F :=
if s is x :: l return 'M_((size_sum s).+1)
then diag_block_mx_rec x l F else 0.

where  is a sequence of natural numbers and is a sequence of algebraic multiplicity of the eigenvalues of a matrix. The function
defines a block in the block diagonal matrix. It takes the size of a block and the location of the block in the block diagonal matrix and returns a matrix of size $n+1$. The block diagonal matrix is defined recursively using the following definition
Fixpoint diag_block_mx_rec k (s : seq nat)
(F : (forall n, nat -> 'M[R]_n.+1)) :=
if s is x :: l return 'M_((size_sum_rec k s).+1)
then block_mx (F k 0
else F k 0

The definition  generates a matrix using four block matrices.  defines the size of a sequence  recursively.

§.§ Complex modulus

We define the modulus of a complex number as
Definition C_mod (x: R[i]):= sqrt ( (Re x)^+2 + (Im x)^+2).

Here, the type  is denoted by , and  and  denote the real and imaginary part of , respectively.

§.§ Matrix and vector norms

In this work, we formalize the 2-norm of a matrix and its induced vector norm. In Coq, we define 2-norm of a matrix as
Definition matrix_norm (n:nat) (A: 'M[complex R]_n.+1) :=
Lub_Rbar (fun x=>
exists v: 'cV[complex R]_n.+1, v != 0 /                x = (vec_norm_C  (mulmx A v))/ (vec_norm_C v)).

where  is the 2-norm of a complex vector, which we define in Coq as
Definition vec_norm_C (n:nat) (x: 'cV[complex R]_n.+1):=
sqrt ([+

The definition  is the least upper bound and is already defined in the  [9] library. Mathematically,  formalizes the following definition of a matrix norm

\begin{equation*}
||A||_i = \sup_{x \neq 0}\frac{||Ax||}{||x||}
\end{equation*}

for a given vector norm $||.||$, which in this case is the L2 vector norm and is mathematically defined as

\begin{equation*}
||x|| = \sqrt{\sum_{j=1}^n |x_j|^2}
\end{equation*}

We will next discuss our formalization of the necessary and sufficient conditions for iterative convergence, which builds on the above mentioned libraries that we discussed.

§ NECESSARY AND SUFFICIENT CONDITIONS FOR ITERATIVE CONVERGENCE

Consider a linear differential equation of the form $\mathcal{A}  u = f$. Here, $\mathcal{A}$ is a linear differential operator acting on a real valued continuous function $u:\mathbb{R}$. For instance, consider a linear differential equation $\frac{d^2u}{dx^2} = 1$. The differential operator $\mathcal{A} = \frac{d^2}{dx^2}$ and $u$ is the real valued function. This differential equation is then discretized in a finite computation domain using a choice of numerical methods. Using a finite difference numerical method to discretize the above-mentioned differential equation in a grid of uniform spacing $h$, we obtain a system of difference equation $A_h u_h = f_h$. Here, the coefficient matrix $A_h$ is known to us and its form depends on the choice of a finite difference scheme. The vector $f_h$ is also known and represents the set of data points. The vector $u_h$, also known as the solution vector is unknown. We therefore solve the linear system $A_h u_h = f_h$ for the solution vector $u_h$. A naive way to solve linear system of equations is to invert the matrix $A_h$ and obtain the solution vector $u_h$ defined as $u_h \overset{\Delta}{=}A_h^{-1}f_h$. However, the matrix $A_h$ could be full and large. The computational complexity of inverting a matrix directly is  $\mathcal{O}(N^3)$, where $N$ is the dimension of the square matrix $A_h$. On the contrary, classical iterative methods like the Jacobi [5] and the Gauss–Seidel iteration [5] have computational complexity of $\mathcal{O}(N^2)$. Therefore, for large systems, iterative methods save significant computational time and are thus commonly used in scientific computing solvers.

From this point forward, we will denote the matrix $A_h$ as $A$, the solution vector $u_h$ as $x$, and the data vector $f_h$ as $b$. The following presentation is a summary of the general theory of iterative methods described in [5].

Let $x$ be a direct solution obtained by inverting the linear system $Ax = b$ as

\begin{equation}\label{direct_solution}
x \overset{\Delta}{=} A^{-1}b
\end{equation}

Here, the matrix $A$ and the vector $b$ are known to us and we are computing the unknown vector $x$. We assume that the matrix $A$ is non-singular. Thus, there exist a unique solution $x$ of the linear system $Ax = b$. For any iterative algorithm, we start with an initial guess vector $x_o$ and obtain a sequence of numerical solutions which are an approximation of the solution $x$. Let $x_m$ be the iterative solution obtained after $m$ iterations obtained by solving the iterative system

\begin{equation}\label{iterative_system}
A_1 x_m + A_2 x_{m-1} = b
\end{equation}

for some choice of initial solution vector $x_o$. The vector $x_{m-1}$ is the iterative solution obtained after $m-1$ iterations. At the $m^{th}$ iteration step, $x_{m-1}$ is known to us. The matrices $A_1$ and $A_2$ are obtained by splitting the original coefficient matrix $A$ such that $A_1$ is easily invertible. Therefore,

\begin{equation}\label{matrix_split}
A = A_1 + A_2
\end{equation}

The choice of matrices $A_1$ and $A_2$ define the choice of an iterative method. For instance, if we choose $A_1$ to be the diagonal entries of matrix $A$ and $A_2$ to be the strictly lower and upper triangular entries of $A$, we obtain the Jacobi method. We discuss about the Jacobi method in detail in Section <ref>. If we choose $A_1$ to be the lower triangular entries of $A$ and $A_2$ to be the strictly upper triangular entries of $A$, we get the Gauss–Seidel iterative method. We discuss about the Gauss–Seidel method in detail in Section <ref>. Therefore, the matrices $A_1$ and $A_2$ are also known to us based on the choice of an iterative method. The data vector $b$ is also known to us. Thus, the unknown solution vector $x_m$ at $m^{th}$ step is obtained by rearranging terms in the iterative system (<ref>) as

\begin{equation}\label{iterative_solution}
x_m = (-A_1^{-1}A_2)x_{m-1} + A_1^{-1}b
\end{equation}

The quantity $-A_1^{-1}A_2$ in equation (<ref>) is called an iterative matrix and we will denote it as $S$. Therefore,

\begin{equation}\label{iterative_matrix}
S \overset{\Delta}{=} -A_1^{-1}A_2
\end{equation}

The iterative convergence error after $m$ iterations is defined as

\begin{equation}
e^m_{iterative} \overset{\Delta}{=}
x_m - x
\end{equation}

The iterative solution $x_m$ is said to converge to $x$ if and only if

\begin{equation}
\lim_{m \to \infty} ||e^m_{iterative}|| = \lim_{m \to \infty} ||x_m - x|| = 0
\end{equation}

where $||.||$ denotes a vector norm. In this work, we will be using the L2 vector norm defines as

\begin{equation*}
||x||_2 = \sqrt{\sum_{j=1}^n |x_i|^2}.
\end{equation*}

The following theorem provides necessary and sufficient conditions for iterative convergence

Let an iterative matrix be defined as (<ref>) for the iterative system (<ref>). The sequence of iterative solutions $\{x_m\}$ converges to the direct solution $x$ for all initial values $x_o$, if and only if the spectral radius of the iterative matrix $S$ is less than $1$ in magnitude.

Spectral radius of a matrix is defined as the maximum of the magnitude of its eigenvalues. We next discuss the proof of Theorem <ref> followed by its formalization in the Coq proof assistant. It is noteworthy that while such proofs have been discussed in numerical analysis literature, we found several missing pieces during the formalization. Most facts about intermediate steps in the proof have just been stated in the numerical analysis literature without a rigorous proof. We therefore spent considerable time developing those proofs during our formalization. In that regard, a contribution of this work is to provide a clean machine-checked proof of the main theorem and any intermediate lemma or fact that was required to close the proof of the Theorem <ref>.

§.§ Proof of Theorem <ref>

To prove Theorem <ref>, we first need to obtain a recurrence relation for the iterative convergence error at $m^{th}$ step in terms of the initial iteration error $(x_o - x)$. Therefore,

\begin{align}
x_{m}-x &= -A_{1}^{-1}A_2 x_{m-1}+A_{1}^{-1}b -x \nonumber \\
&= -A_{1}^{-1}A_{2}x_{m-1}+ A_{1}^{-1}(Ax) -x \quad [Ax\overset{\Delta}{=}b] \nonumber \\
&= -A_{1}^{-1}A_{2}x_{m-1}+A_{1}^{-1}(A_{1}+A_{2})x -x \qquad [A_{1}+A_2 = A] \nonumber \\
&=-A_{1}^{-1}A_{2}x_{m-1}+A_{1}^{-1}A_1x +A_{1}^{-1}A_{2}x -x \nonumber \\
&=-A_{1}^{-1}A_{2}(x_{m-1}-x) \quad [A_{1}^{-1}A_{1}\overset{\Delta}{=}I]\nonumber
\end{align}

Taking norm of the vector on both sides, the iterative convergence error at the $m^{th}$ step can be written in terms of the iterative convergence error at $(m-1)^{th}$ step as

\begin{equation}\label{iter1}
||x_m - x|| = || (-A_1^{-1}A_2) (x_{m-1} -x) ||
\end{equation}

Since, the system is linear, equation (<ref>) can be written in terms of the initial iteration error as

\begin{equation}\label{iter2}
||x_m - x|| = || (-A_1^{-1}A_2)^m (x_o - x)||
\end{equation}

Taking limits of the vector norms on both sides of equation (<ref>),

\begin{equation}\label{iter3}
\lim_{m \to \infty} ||x_m - x|| = \lim_{m \to \infty} ||(-A_1^{-1}A_2)^m (x_o - x)||
\end{equation}

If $x_o = x$, the iterative convergence error is zero trivially. The case $x_o \neq x$ is interesting and we can prove Theorem <ref> by splitting it into two parts

\begin{equation}\label{subproof_1}
(\forall x_o,  \lim_{m \to \infty} ||(-A_1^{-1}A_2)^m (x_o - x)|| = 0) \iff \lim_{m \to \infty} ||(-A_1^{-1}A_2)^m|| = 0
\end{equation}

\begin{equation}\label{subproof_2}
\lim_{m \to \infty} ||(-A_1^{-1}A_2)^m|| = 0 \iff \rho(-A_1^{-1}A_2) < 1
\end{equation}

where $||(-A_1^{-1}A_2)^m||$ is the matrix norm. We will be considering 2-norm of matrix and its induced vector norm in our work. The quantity $\rho(-A_1^{-1}A_2)$ is the spectral radius of the iteration matrix, $S = (-A_1^{-1}A_2)$ and is defined as

\begin{equation*}
\rho(S) = \max_{i} \{|\lambda_i(S)|\}, \forall i = 0 .. N-1
\end{equation*}

where $\lambda(S)_i$ is the $i^{th}$ eigenvalue of $S$. Therefore,

\begin{equation*}
\rho(S) < 1 \iff (\forall i, i < N \implies |\lambda_i(S)|)
\end{equation*}

We next discuss the proofs of (<ref>) and (<ref>).

§.§.§ Proof of (<ref>)

We need to prove that

\begin{equation*}
\lim_{m \to \infty} ||(-A_{1}^{-1}A_{2})^m|| = 0 \implies  (\forall x_o,  \lim_{m \to \infty} ||(-A_1^{-1}A_2)^m (x_o - x)|| = 0)
\end{equation*}

: Given $x_o$,

\begin{align}
\lim_{m \to \infty} ||(-A_{1}^{-1}A_{2})^{m} (x_{o}-x)||
&\leq \lim_{m \to \infty}||(-A_{1}^{-1}A_{2})^m||\,||x_{o}-x||;\quad [||Ax||\leq ||A||||x||]\nonumber \\
&= \left(\lim_{m \to \infty}||(-A_{1}^{-1}A_{2})^m|| \right) \left(\lim_{m \to \infty}||x_{o}-x|| \right) \nonumber \\
&=0;  \quad [\text{since, } \lim_{m \to \infty}||(-A_{1}^{-1}A_{2})^m||=0]\nonumber
\end{align}

We need to prove that

\begin{equation*}
(\forall x_o,  \lim_{m \to \infty} ||(-A_1^{-1}A_2)^m (x_o - x)|| = 0) \implies \lim_{m \to \infty}||(-A_{1}^{-1}A_{2})^m||=0
\end{equation*}

We start by unfolding the definition of the norm of the iterative matrix

\begin{equation}\label{matrix_norm}
||(-A_1^{-1}A_2)^m|| = \sup_{(x_o - x) \neq 0} \frac{||(-A_1^{-1}A_2)^m (x_o - x) ||}{||x_o - x||}
\end{equation}

Therefore, we need to prove that

\begin{equation}\label{limit_matrix_norm}
\lim_{m \to \infty}
\left(\sup_{(x_o - x) \neq 0}
\frac{||(-A_1^{-1}A_2)^m (x_o - x) ||}{||x_o - x||} \right) = 0
\end{equation}

We can prove (<ref>) by choosing an upper bound for the matrix norm in (<ref>), proving that the limit of this upper bound converges to zero and then applying the sandwich theorem [26] for limits. We choose this upper bound as $\sum_{j < N} ||(-A_1^{-1}A_2)^m e_j||$, i.e.,

\begin{equation}\label{upper_bound}
\sup_{(x_o - x) \neq 0} \frac{||(-A_1^{-1}A_2)^m (x_o - x) ||}{||x_o - x||} \leq \sum_{j < N} ||(-A_1^{-1}A_2)^m e_j||
\end{equation}

where $e_j$ is the unit vector corresponding to a principal direction in the cartesian coordinate system, i.e.,
$e_j = \textbf{1}_j$. The vector $\textbf{1}_j$ is a unit vector with the entry $1$ in the $j^{th}$ place and other entries in the vector being $0$.

\begin{align*}
&\lim_{m \to \infty}
\left(\sup_{(x_o - x) \neq 0}
\frac{||(-A_1^{-1}A_2)^m (x_o - x) ||}{||x_o - x||} \right) \leq
\lim_{m \to \infty}
\sum_{j < N} ||(-A_1^{-1}A_2)^m e_j|| \\
\lim_{m \to \infty}
\left(\sup_{(x_o - x) \neq 0}
\frac{||(-A_1^{-1}A_2)^m (x_o - x) ||}{||x_o - x||} \right)
\leq \sum_{j < N}
\left(\lim_{m \to \infty}
||(-A_1^{-1}A_2)^m e_j||\right)
\end{align*}

We can quantify $x_o$ in the hypothesis with $x + e_j, \forall j, j < N$. Therefore, $\forall j, j< N$, we have from the hypothesis,

\begin{equation*}
\lim_{m \to \infty}
|| (-A_1^{-1}A_2) e_j|| = 0
\end{equation*}

\begin{equation*}
\sum_{j < N}
\left(\lim_{m \to \infty}
||(-A_1^{-1}A_2)^m e_j||\right) = 0
\end{equation*}

We can then apply the sandwich theorem [26] for limits to prove that

\begin{equation}\label{limit_matrix_norm_0}
\lim_{m \to \infty}
\left(\sup_{(x_o - x) \neq 0}
\frac{||(-A_1^{-1}A_2)^m (x_o - x) ||}{||x_o - x||} \right) = 0
\end{equation}

We justify the choice of the upper bound in (<ref>) in the following proof.
We can decompose the vector $x_o - x$ into its components along the principal axes in the cartesian coordinate system as

\begin{equation*}
x_o - x = \sum_{j < N} (x_o - x)_j e_j
\end{equation*}

\begin{equation}\label{linear_compose}
(-A_1^{-1}A_2)^m (x_o - x) = \sum_{j < N} (-A_1^{-1}A_2)^m (x_o - x)_j e_j
\end{equation}

By taking a vector norm on both sides of (<ref>),

\begin{equation*}
||(-A_1^{-1}A_2)^m (x_o - x)|| \leq
||\sum_{j < N} (-A_1^{-1}A_2)^m (x_o - x)_j e_j||
\end{equation*}

Using the triangle inequality property of the vector norm, we get

\begin{equation}\label{ineq_norm}
||(-A_1^{-1}A_2)^m (x_o - x)|| \leq
\sum_{j < N} ||(-A_1^{-1}A_2)^m (x_o - x)_j e_j||
\end{equation}

Since $(x_o - x) \neq 0 $, $||x_o - x|| \neq 0$. Hence, dividing by $||x_o - x||$ on both sides of (<ref>),

\begin{align*}
\frac{||(-A_1^{-1}A_2)^m (x_o - x)||}{||x_o - x||} &\leq
\sum_{j < N} \frac{||(-A_1^{-1}A_2)^m (x_o - x)_j e_j||} {||x_o -x||}\\
\sum_{j < N}
\frac{|(x_o - x)_j|}{||x_o -x||}
||(-A_1^{-1}A_2)^m e_j|| \\
&\leq \sum_{j < N}
||(-A_1^{-1}A_2)^m e_j||
\end{align*}

Here, we first use the absolute homogenity property ($||ax|| = |a| ||x||$, for any scalar $a$ and vector $x$) of the vector norm. Then we use the fact that

\begin{equation*}
|(x_o - x)_j| \leq
||x_o - x||
\end{equation*}

§.§.§ Proof of (<ref>)

(Sufficiency): We need to prove that

\begin{equation*}
\lim_{m \to \infty} ||(-A_{1}^{-1}A_2)^m||=0 \implies \rho(-A_{1}^{-1}A_2) < 1
\end{equation*}

\begin{equation*}
\rho(-A_{1}^{-1}A_2) =\max_{0 \leq i < n}|\lambda_i (-A_{1}^{-1}A_2)|,
\end{equation*}

\begin{equation*}
\lim_{m \to \infty} ||(-A_{1}^{-1}A_2)^m||=0 \implies
(\forall i, 0 \leq i < n, \; |\lambda_i|<1)
\end{equation*}

\begin{equation*}
\lim_{m \to \infty} |\lambda_{i}|^m =0 \implies |\lambda_i|<1, \; \forall i, 0 \leq i < n
\end{equation*}

to the goal statement, we need to prove:

\begin{equation*}
\lim_{m \to \infty} |\lambda_i|^m = 0, \;\forall i, \; 0 \leq i < n
\end{equation*}

under the hypothesis $\lim_{m \to \infty} ||(-A_{1}^{-1}A_2)^m||=0$. We now use the definition of an eigensystem, i.e.,

\begin{equation*}
(-A_{1}^{-1}A_2) v_{i} = \lambda_i v_i
\end{equation*}

where $v_i$ is an eigenvector corresponding to the eigenvalue $\lambda_i$. Therefore,

\begin{align}
\lim_{m \to \infty} |\lambda_i|^m = 0 &\implies
\lim_{m \to \infty} |\lambda_i|^m ||v_i|| = 0 \nonumber \\
& \implies \lim_{m \to \infty} ||\lambda_i ^m v_i||=0;
\quad [|\lambda_i|^m = |\lambda_i^m| \; \land\; |\lambda_i^m| ||v_i|| = ||\lambda_i ^m v_i||]
\nonumber \\
& \implies \lim_{m \to \infty} ||(-A_{1}^{-1}A_2)^m v_{i}||=0; \quad [(-A_{1}^{-1}A_2)^m v_{i} = \lambda_i ^m v_i] \label{part2_suf_1}
\end{align}

But the compatibility relation for vector and matrix norms dictates,

\begin{equation}\label{part2_suf_2}
0 \leq ||(-A_{1}^{-1}A_2)^m v_{i}|| \leq ||(-A_{1}^{-1}A_2)^m|| ||v_i||
\end{equation}

Since, $||v_i|| \neq 0 $ by definition of an eigensystem,

\begin{equation}
\lim_{m \to \infty} ||(-A_{1}^{-1}A_2)^m||= 0 \implies \lim_{m \to \infty} ||(-A_{1}^{-1}A_2)^m|| ||v_i|| =0
\end{equation}

We can then apply the sandwich theorem [26] to prove that

\begin{equation*}
\lim_{m \to \infty} || (-A_1^{-1}A_2)^m v_i|| = 0
\end{equation*}

We need to prove that

\begin{equation*}
\rho(-A_{1}^{-1}A_2) <1 \implies \lim_{m \to \infty} ||(-A_{1}^{-1}A_2)^m||=0.
\end{equation*}

To prove necessity, we obtain the Jordan decomposition of the iterative matrix, $S=(-A_{1}^{-1}A_2)$. From the Jordan normal form theorem [5], we know that there exists $V, J \in \mathbb{C}^{n \times n}$, $V$ non-singular and $J$ block diagonal such that:

\begin{equation}\label{jordan_form}
S= V J V^{-1}
\end{equation}

\begin{equation*}
J =
\begin{bmatrix}
J_{k_1}(\lambda_1) & 0 & 0 & \hdots & 0\\
0 & J_{k_2}(\lambda_2) & 0 & \hdots & 0 \\
\vdots & \hdots & \ddots & \hdots & \vdots \\
0 & \hdots & 0 & J_{k_{s-1}}(\lambda_{s-1}) & 0\\
0 & \hdots & \hdots & 0 & J_{k_s}(\lambda_{s})
\end{bmatrix};
\quad
J_{k_i}(\lambda_i) =
\begin{bmatrix}
\lambda_i & 1 & 0 & \hdots & 0\\
0 & \lambda_i & 1 & \hdots & 0 \\
\vdots & \vdots & \ddots & \ddots & \vdots \\
0 & 0 & \hdots & \lambda_i & 1\\
0 & 0 & \hdots & 0 & \lambda_i
\end{bmatrix}
\end{equation*}

where $J_{k_i}(\lambda_i)$ is the Jordan block corresponding to the eigenvalue $\lambda_i$.
Thus, $S^m = V J^m V^{-1}$. Since $J$ is block diagonal,

\begin{equation*}
J^m =
\begin{bmatrix}
J_{k_1}^m(\lambda_1) & 0 & 0 & \hdots & 0\\
0 & J_{k_2}^m(\lambda_2) & 0 & \hdots & 0 \\
\vdots & \hdots & \ddots & \hdots & \vdots \\
0 & \hdots & 0 & J_{k_{s-1}}^m(\lambda_{s-1}) & 0\\
0 & \hdots & \hdots & 0 & J_{k_s}^m(\lambda_{s})
\end{bmatrix}
\end{equation*}

Now, a standard result on the $m^{th}$ power of a $k_i \times k_i$ Jordan block states that, for $m \geq k_i-1$

\begin{equation*}
J_{k_i}^m(\lambda_i) =
\begin{bmatrix}
\lambda_i^m & {m \choose 1}\lambda_i^{m-1} & {m \choose 2}\lambda_i^{m-2} & \hdots & {m \choose {k_i-1}}\lambda_i ^{m - k_i +1} \\
0 &  \lambda_i^m & {m \choose 1}\lambda_i^{m-1} & \hdots &
{m \choose {k_i-2}}\lambda_i ^{m - k_i +2}\\
\vdots & \vdots & \ddots & \ddots & \vdots\\
0 & 0 & \hdots & \lambda_i^m & {m \choose 1}\lambda_i^{m-1} \\
0 & 0 & \hdots & 0 & \lambda_i^m
\end{bmatrix}
\end{equation*}

To prove that $\lim_{m \to \infty} ||S^m||=0$, we need to prove that $\lim_{m \to \infty} ||J^m|| = 0 $, since $V$ is non-singular, i.e. $||V|| \neq 0$. Since the 2-norm of a matrix is bounded above by the Frobenius matrix norm as:
$ ||J^m||_2 \leq ||J^m||_F
we can prove $\lim_{m \to \infty} ||J^m||_2 =0 $ by proving  $\lim_{m \to \infty} ||J^m||_F = 0$. The Frobenius norm of matrix $J^m$ is defined as:

\begin{equation*}
||J^m||_F \overset{\Delta}{=} \sqrt{\sum_{i} \sum_{j} |(J^m)_{(i,j)}|^2}
\end{equation*}

Thus, we need to prove that

\begin{equation}\label{each_entry_zero}
\forall i \; j, \lim_{m \to \infty} |(J^m)_{(i,j)}|^2 = 0
\end{equation}

The zero entries of the Jordan block diagonal matrix tend to zero trivially.
Since, $|\lambda_i| < 1$, the diagonal elements of the Jordan block, $J_{k_i}^m (\lambda_i)$ tend to zero. Proving that the off diagonal elements tend to zero is more difficult.

We have to prove that $\lim_{m \to \infty} {m \choose k} |\lambda_i |^{m-k} = 0, \; \forall k \leq k_i -1$. The function $f(m) = {m \choose k}$ increases as $m$ increases while the function $g(m) = |\lambda_i|^{m-k}$ decreases as $m$ increases since $0 < |\lambda_i| < 1$. Informally, it can be argued that $g(m)$ decreases at a faster rate than the function $f(m)$. Hence, the limit should tend to zero. But proving it formally in Coq was challenging. We need to first bound the function $f(m)$ with a function $h(m) = \frac{m^k} {k!}$ to obtain a sequence, $\frac{m^k}{k!}|\lambda_i|^{m-k}$, for which it would be easy to prove the limit. Therefore, we split the proof into proofs of two facts:

\begin{equation}\label{choose_prove}
{m \choose k} \leq \frac{m^k}{k!}
\end{equation}

\begin{equation}\label{sequence_prove}
\lim_{m \to \infty} \frac{m^k}{k!}|\lambda_i|^{m-k}=0
\end{equation}

An informal proof of the inequality (<ref>) is available in Appendix <ref>.

In order to prove $\lim_{m \to \infty} \frac{m^k}{k!} |\lambda_i|^{m-k} = 0$, we use the ratio test for convergence of sequences. The formalization of ratio test has not yet been done in Coq to our knowledge. So, we formalized the ratio test since it provides an easier test for proving convergence of sequences as compared to first bounding the function with an easier function for which the convergence could be proved with the existing Coq libraries. The process of bounding the function $\Gamma (m) = f(m)g(m)$ with an easier function to prove convergence was challenging for us due to the behavior of $\Gamma (m)$. Using plotting tools like Wolfram plot or MATLAB plot, we observed that for

$|\lambda_i| \leq 0.5$, the function $\Gamma (m)$ was monotonously decreasing, while for $0.5 < |\lambda_i| < 1 $, $\Gamma (m)$ increases first and then decreases. Moreover, the location of the maxima of $\Gamma (m) $ in the interval $0.5 < |\lambda_i| < 1$ depends of the number of iterations. Hence, bounding $\Gamma (m)$ with a monotonously decreasing function is challenging for $0.5 < |\lambda_i| < 1$. But we know that $\Gamma (m) $ decays eventually. For such scenarios, just comparing the terms $\Gamma (m+1) $ and $\Gamma (m)$ provides a much simpler test for the convergence of the sequence $\Gamma (m)$. This is where ratio test for the convergence of sequence comes to rescue. The ratio test for the convergence of sequences is stated as follows.

[27]
If $(a_n)$ is a sequence of positive real numbers such that $\lim_{n \to \infty} \frac{a_{n+1}}{a_n} = L$ and $L < 1$, then $(a_n)$ converges and $\lim_{n \to \infty} a_n = 0$.

We provide an informal proof of Theorem <ref> in the Appendix <ref>.
In our case, the sequence $a_n = n^k |\lambda_i|^n$. Therefore, the ratio $\frac{a_{n+1}}{a_n} = \frac{(n+1)^k|\lambda_i|^{n+1}}{n^k |\lambda_i|^n}$. Therefore,

\begin{align*}
\lim_{n \to \infty} \frac{a_{n+1}}{a_n} &= \frac{(n+1)^k |\lambda_i|^{n+1}}{n^k |\lambda_i|^n}\\
&= \lim_{n \to \infty}\left( 1+ \frac{1}{n}\right)^k |\lambda_i|  \\
& = |\lambda_i| \left[\lim_{n \to \infty}\left( 1+ \frac{1}{n}\right)^k\right] \\
&= |\lambda_i| \left[\lim_{n \to \infty}\left( 1+ \frac{1}{n}\right)\right]^k  \\
\end{align*}

Therefore, $L=|\lambda_i|$, which we know is less than $1$. Thus, we can apply Theorem <ref> to prove that $\lim_{n \to \infty} n^k |\lambda_i|^n = 0$. Since $k!$ and $|\lambda_i|^k$ are non zero constants,
$ lim_{m \to \infty} \frac{m^k}{k!} |\lambda_i|^{m-k} = 0. $

§.§ Formalization in Coq

We state Theorem <ref> in Coq as follows:
Theorem iter_convergence:
forall (n:nat) (A: 'M[R]_n.+1) (b: 'cV[R]_n.+1) (A1 A2 : 'M[R]_n.+1),
A ∈unitmx ->
A1 ∈unitmx ->
A = A1 + A2 ->
let x := (invmx A) *m b in
(let S_mat:= RtoC_mat (- ( A1^-1 *m A2)) in
(forall (i: 'I_n.+1), (C_mod (lambda S_mat i) < 1)
(forall x0: 'cV[R]_n.+1,
is_lim_seq (fun m:nat => vec_norm ((X_m m.+1 x0 b A1 A2) - x)) 0

Since we deal with a generic case where a real matrix is allowed to have complex eigenvalues and eigenvectors, we also work in the complex field.  transforms a real matrix to a complex matrix so as to be consistent with types. In , a complex number is defined on top of a real closed field [28, 29]. Thus, given a real matrix $A$,  transforms a real entry $A_{ij}: \mathbb{R}$ to a complex number $\tilde{A}_{ij}:\mathbb{C}:= (A_{ij} +i* 0)\%C$. In Coq, we formally define  as follows:
Definition RtoC_mat (n:nat) (A: 'M[R]_n):'M[complex R]_n :=
_(i<n, j<n) ((A i j) +i* 0)

Let us now discuss the hypothesis of the theorem statement. The first hypothesis states that the matrix $A$ is invertible. Hence, there exists a unique solution to the linear system $Ax= b$. We define this solution vector $x$ in Coq using the let binding
let x := (invmx A) *m b

where the  definition  defines the inverse of a matrix. The hypothesis  corresponds to the definition (<ref>). As discussed earlier, we want the matrix $A_1$ to be easily invertible so as to construct the iterative matrix $S= -(A_1^{-1}A_2)$. The invertibility condition is stated by the hypothesis
A1 ∈unitmx

We define the iterative solution after $m$ steps, $x_m$ from the iterative system (<ref>) as
Fixpoint X_m (m n:nat) (x0 b: 'cV[R]_n.+1) (A1 A2: 'M[R]_n.+1) : 'cV[R]_n.+1:=
match m with
| O => x0
| S p => ((- ((A1^-1) *m A2)) *m (X_m p x0 b A1 A2)) + ((A1^-1) *m b)

We define the iterative matrix $S_{mat}$ in the let binding of the theorem statement in Coq.

In our formalization of the iterative convergence, we rely on the existing formalization of the Jordan canonical forms by Guillaume Cano and Maxime  Dénès [13]. We use their definition of eigenvalues of a matrix derived from the roots of the characteristic polynomials of the Smith Normal form of a matrix $A$.

We first define a sequence of eigenvalues as the diagonal entries of Jordan matrix.
[language = Coq]
Definition lambda_seq (n: nat) (A: 'M[complex R]_n.+1):=
let sizes:= size_sum
[seq x.2.-1
| x <- root_seq_poly (invariant_factors A)] in
[seq (Jordan_form A) i i | i <- enum 'I_sizes.+1].

returns a sequence of pair of roots and its multiplicity, of the polynomial . The  are the polynomials in the diagonal of the Smith Normal form of a matrix. In this case, the sequence contains the pair of eigenvalues of matrix A and its multiplicity. The Jordan form matrix contains these eigenvalues in its diagonal which we extract using .

The $i^{th}$ eigenvalue of a matrix $A$ is then defined as the $i^{th}$ component of the sequence of eigenvalues .
[language = Coq]
Definition lambda (n: nat) (A: 'M[complex R]_n.+1) (i: 'I_n.+1) :=
(@nth _  0

To take full advantage of the lemmas describing eigenvalues and eigenvectors as defined in , we had to relate the definition of eigenvalue, , and the one defined in . In , an eigenvalue $a$, of a matrix $A$ is defined as a predicate
[language = Coq]
Definition eigenvalue : pred F := fun a => eigenspace a != 0.

stating that the eigenspace corresponding to this eigenvalue is non-zero.
[language = Coq]
Lemma Jordan_ii_is_eigen:
forall (n: nat) (A: 'M[complex R]_n.+1),
forall  (i: 'I_n.+1),
@eigenvalue (complex_fieldType _) n.+1 A (@nth _  0

The lemma  asserts that the  that we obtain from the Jordan matrix of $A$ is an eigenvalue of $A$ using the  predicate in .

It should be noted that a square matrix $A$ is assumed to be of size at least 1. This design choice is justified by Cano and Dénès in [13]. A square matrix only forms a ring when their size is convertible to the successor of an integer [13]. The ring property of a square matrix was really helpful in the formalization of diagonal block matrices by Cano and Dénès. Therefore, to use their formalization of canonical forms, we also stick with denoting the type of a complex matrix as . Here, we instantiate a generic ring type with a complex ring type.

As discussed earlier in the section on informal proof, we split the proof of the Theorem <ref> into two sub proofs. The statement of these sub proofs are formalized in Coq as:

(Part 1):
(forall x0: 'cV[R]_n.+1,
is_lim_seq (fun m : nat => vec_norm  ((X_m m.+1 x0 b A1 A2) - x)) 0
is_lim_seq (fun m:nat =>
(matrix_norm  (RtoC_mat ((- (A1^-1 *m A2))^+m.+1) ))) 0

(Part 2):
is_lim_seq (fun m:nat => matrix_norm
(RtoC_mat ((- ((A1^-1 *m A2)))^+m.+1 ))) 0
(let S_mat := RtoC_mat (- (A1^-1 *m A2)) in
forall i : 'I_n.+1, (C_mod (lambda S_mat i) < 1)
We have used the  predicate from the  [9] library to define the limits for the sequence of solution vectors $x_{m}$ converging to $x$.

For the proof of the fact in the forward direction in Part 1, we prove the following lemma in Coq
Lemma lim_max: forall (n:nat) (A: 'M[R]_n.+1) (x: 'cV[R]_n.+1),
(forall x0: 'cV[R]_n.+1,
let v:= x0 - x in
let vc:= RtoC_vec v in
is_lim_seq (fun m: nat => vec_norm_C ((RtoC_mat (A^+m.+1)) *m vc)) 0
is_lim_seq (fun m:nat => matrix_norm (RtoC_mat (A^+m.+1))) 0

The lemma  corresponds to the proof of (<ref>). As discussed earlier, to prove , an important step was to decompose the vector $x_o - x$ into its components along the principal axis of the Cartesian coordinate system. We prove this fact in Coq using the following lemma statement
Lemma base_exp: forall (n:nat) (x: 'cV[complex R]_n.+1),
x = [+

where we defined the principal unit vector  in Coq as
Definition e_i n:nat (i: 'I_n.+1): 'cV[complex R]_n.+1 :=
_(j < n.+1) (if (i= =j :> nat) then (1 +i* 0)

To prove the sufficiency in part 2, we had to prove the following lemma:
Lemma is_lim_seq_geom_nec (q:R):
is_lim_seq (fun n => (q ^ n.+1)

While a lemma exists for other direction in the  [9] formalization of limits, lemma  was missing.

As discussed earlier, to prove sufficiency condition for iterative convergence, we had to formalize the ratio test for convergence of sequences which was missing in the existing Coq libraries.
In Coq, we state Theorem <ref> as:
Lemma ratio_test: forall (a: nat -> R) (L:R),
(0 < L / L < 1) ->
(forall n:nat, (0 < a n)
(is_lim_seq ( fun n:nat =>  ((a (n.+1))/(a n))
is_lim_seq (fun n: nat => a n) 0

We then use the lemma  to formally prove

\begin{equation*}
\lim_{m \to \infty} (m+1)^k x^{m+1} = 0
\end{equation*}

which we state in Coq as:
Lemma lim_npowk_mul_to_zero: forall (x:R) (k:nat),
(0 < x)
is_lim_seq (fun m:nat => ((m.+1)

To bound the binomial coefficient ${m \choose k}$ with $\frac{m^k}{k!}$, we formally state the inequality (<ref>) as follows
Lemma choice_to_ring_le: forall (n k:nat),  (0<n)

Since the inequality holds in a complex field, we have to first convert the complex inequality to a real inequality. This was possible using the existing formalization of complex fields in . Since the choice function and the factorials are naturals in Coq, we then define a lemma which proves that if the inequality holds for naturals, the inequality also holds in real field after a coercion is defined from naturals to reals. In Coq, we formalize the lemma as:
[language= Coq]
Lemma nat_ring_mn_le: forall (m n:nat), (m<= n)

Applying the lemma , we obtain the following inequality

\begin{equation}\label{choice_nat}
\left({n \choose k} \leq \frac{n^k} { k!}\right)\% N
\end{equation}

It should be noted that the division in (<ref>) is an integer division. The inequality (<ref>) can easily be proved as discussed in the previous section. However, to use the inequality (<ref>) in the proof of lemma , one needs to be careful with converting an integer division to a real division.

To prove that each element of the Jordan block matrix converges to zero as in equation (<ref>), we prove the following lemma in Coq
Lemma each_entry_zero_lim:
forall (n:nat) (A: 'M[complex R]_n.+1),
let sp := root_seq_poly (invariant_factors A) in
let sizes := [seq x0.2.-1 | x0 <- sp] in
forall i j: 'I_(size_sum sizes).+1,
(forall i: 'I_(size_sum sizes).+1 ,
(C_mod (@nth _ (0 +i* 0)
(fun m: nat =>
(C_mod (diag_block_mx sizes
(fun n0 i1 : nat =>
_(i2, j0)
(('C(m.+1, j0 - i2))
(nth (0, 0
(m.+1 - (j0 - i2)) *+ (i2 <= j0))) i j)) ²) 0

The lemma  states that if the magnitude of each eigenvalue of a matrix $A$ is less than $1$, i.e., $|\lambda_i(A)| < 1, \; \forall i, \; 0 \leq i < N$, then the limit of each term in the expanded Jordan matrix is zero as $m$ approaches $\infty$. Here, the block diagonal matrix  takes an expanded Jordan block $J_{k_i}^m(\lambda_i), \forall i, 0\leq i < N$ and constructs the Jordan matrix $J^m$. We then take a modulus of each entry of $J^m$ and prove its limit being zero as $m \to \infty$.

A key challenge we faced when proving the above lemma was extracting each Jordan block of the diagonal block matrix. The diagonal block matrix is defined recursively over a function which takes a block matrix of size $\mu_i$ denoting the algebraic multiplicity of each eigenvalues $\lambda_i$. We had to carefully destruct this definition of diagonal block matrix and extract the Jordan block and the zeros on the off-diagonal entries. We can then prove the limit on this Jordan block by exploiting its upper triangular structure.

We state the lemma to extract a Jordan block from the block diagonal matrix as follows:
[language = Coq]
Lemma diag_destruct (R: ringType)
(s : seq nat) (F : (forall n, nat -> 'M[R]_n.+1)):
forall i j: 'I_(size_sum s).+1,
(exists k l m n,
(k <= size_sum s)
(forall p:nat, (diag_block_mx s (fun k l:nat => (F k l)^+ p.+1)) i j =
((F k l)^+ p.+1) m n) /    (diag_block_mx s F i i = (F k l) m m))
(forall p:nat, (diag_block_mx s (fun k l:nat => (F k l)^+ p.+1)) i j = 0).

where  is the sum of the algebraic multiplicities of the eigenvalues and equals the total size of the matrix $n$. We prove this fact using the following lemma statement in Coq
[language = Coq]
Lemma total_eigen_val: forall (n:nat) (A: 'M[complex R]_n.+1),
(size_sum [seq x.2.-1 | x <- root_seq_poly (invariant_factors A)]).+1 = n.+1.

The lemma  helps us get around the dimension constraint imposed by the design of the Jordan form of a matrix $A$.
Limits of the off-diagonal elements can then be trivially proven to zero. This completes the proof of sufficiency condition for convergence of iterative convergence error.

We will next discuss an example problem that we instantiate our iterative convergence theorem with.

§ EXAMPLE PROBLEM

In our work, the linear differential equation $\mathcal{A}u = f$ that we chose was $\frac{d^2 u}{dx^2} = 1$ for $x \in (0,1)$ and the boundary conditions being $u(0) = u(1) = 0$. Here, the differential operator $\mathcal{A}$ is $\frac{d^2}{dx^2}$, and $f$ is the constant function $1$. We chose a uniform grid with $N$ points in the interior of the 1D domain. The grid has a uniform spacing $h$. We will be using a central difference scheme [19] for discretizing the differential equation. Therefore, the difference equation at point $x_i$ in the interior of the 1D domain is given by

\begin{equation}\label{diff}
\frac{-u(x_{i+1})+2u(x_i)-u(x_{i-1})}{h^2}=-1;
\quad h=x_{i+1}-x_i=\frac{1}{N+1}
\end{equation}

When we stack the equation (<ref>) for all points in the interior of the 1D domain, we get a linear matrix system

\begin{equation}\label{matrix_A}
\underbrace{
\frac{1}{h^2}
\begin{bmatrix}
2 & -1 & 0 & 0 & 0 &  \hdots & 0 \\
-1 & 2 & -1 & 0 & 0 & \hdots & 0\\
\vdots & & \ddots & \ddots & \ddots  & & \vdots\\
0 & \hdots & 0 & 0 & -1 & 2 & -1\\
0 & \hdots & 0 & 0 & 0  & -1 & 2
\end{bmatrix}}_{A}
\underbrace{
\begin{bmatrix}
u_1 \\
u_2 \\
\vdots\\
u_N-1 \\
\end{bmatrix}}_{x} =
\underbrace{
\begin{bmatrix}
-1 \\
-1 \\
\vdots \\
\end{bmatrix}}_b
\end{equation}

Here, $A$ is the coefficient matrix, $b$ is the data vector and $x$ is the unknown solution vector, which can be exactly obtained by inverting the matrix $A$, i.e., $x = A^{-1}b$. But, we will obtain an approximation of $x$ using iterative algorithms. In the next couple of sections, we will discuss two classical iterative methods: Gauss–Seidel and Jacobi iteration. We will instantiate these classical iterative algorithms with this example problem and apply Theorem <ref> to prove convergence of the approximate solutions, obtained using these algorithms,  to the exact solution.

§ GAUSS–SEIDEL ITERATION

The Gauss–Seidel iterative method is an instance of iterative methods to solve the system of linear equations
$ Ax = b \quad \text{or} \quad
\sum_{j=1}^{n} A_{ij}x_j = b_i,\; \forall i= 1 .. n
$, using the following update formula

\begin{equation}\label{gauss_system}
x_i^{m+1} = \frac{1}{A_{ii}} \left(
b_i - \sum_{j=1}^{i-1} A_{ij}x_j^{m+1} - \sum_{j=i+1}^{n} A_{ij}x_j^m \right), \quad i = 1, 2, \hdots, n.
\end{equation}

where $x_i^{m+1}$ is the value of $x_i$ after $m+1$ iterations, and $x_i^m$ is the value of $x_i$ after $m$ iterations.

From equation (<ref>), we can construct matrix $A_1$ and $A_2$

as [5]:

\begin{equation}
\begin{cases}
A_{ij} & \text{if $i \geq j$} \\
0 & \text{if $i < j$}
\end{cases};
\qquad
\begin{cases}
A_{ij} & \text{if $i < j$} \\
0 & \text{if $i \geq j$}
\end{cases}
\end{equation}

The following theorem by Reich [30] provides a criteria the convergence of the Gauss–Seidel iterative method.

[30]
If A is real, symmetric nth-order matrix with all terms on its main diagonal positive, then a sufficient condition for all the n characteristic roots of $(-A_{1}^{-1}A_2)$ to be smaller than unity in magnitude is that A is positive definite.

From an application point of view, only the sufficiency condition is important. This is because to apply Theorem (<ref>), we only need to know that the magnitude of the eigenvalues are less than 1.  Thus, to prove the convergence of Gauss–Seidel iteration, we first apply Theorem <ref> to get the eigenvalue condition in the goal and then apply Theorem <ref> to complete the proof. Since computing eigenvalues are not very trivial in most cases, the positive definite property of the matrix $A$ provides an easy test for $|\lambda| < 1$ for Gauss–Seidel iteration matrix. The proof of necessity uses an informal topological argument that would be difficult to formalize.
Therefore, it is enough to just formalize the sufficiency condition.

Next we present an informal proof [30] followed by its formalization in the Coq proof assistant.

Let $z_i$ be the $i^{th}$ characteristic vector of $-(A_{1}^{-1}A_2)$ corresponding to the characteristic root $\mu_{i}$. Then

\begin{equation}
-(A_{1}^{-1}A_{2}) z_i = \mu_{i}z_{i}
\end{equation}

Multiplying by $-(\bar{z_{i}}'A_{1})$ on both sides,

\begin{equation}\label{reich_2}
(-\bar{z_{i}}'A_{1})(-A_{1}^{-1}A_{2})z_i = -\mu_{i}\bar{z_{i}}'A_{1}z_{i}
\end{equation}

where $\bar{z_{i}}'$ is the conjugate transpose of $z_i$ obtained by taking the conjugate of each element of $z_i$ followed by transpose of the vector. Equation (<ref>) then simplifies to:

\begin{equation}
\bar{z_i}'A_2 z_i = -\mu_i \bar{z_i}'A_1 z_i; \quad [A_1 A_{1}^{-1}=I]
\end{equation}

Consider the bi-linear form, $\bar{z_i}'Az_i$,

\begin{equation}
\bar{z_i}'Az_{i}= \bar{z_i}'(A_1+A_2)z_i = \bar{z_{i}}'A_1 z_i + \bar{z_i}'A_2 z_i = \bar{z_{i}}'A_1 z_i - \mu_i \bar{z_i}'A_1 z_i = (1-\mu_i)\bar{z_{i}}'A_1 z_i  \label{reich_3}
\end{equation}

Taking conjugate transpose of equation (<ref>) on both sides,

\begin{equation}\label{reich_4}
\bar{z_i}'Az_{i} = (1-\bar{\mu_i})\bar{z_{i}}'A_1' z_i
\end{equation}

Let $D$ be the diagonal matrix defined as:

\begin{equation}
\begin{cases}
A_{ij} & \text{ if $i=j$}\\
0 & \text{if $i \neq j$}
\end{cases}
\end{equation}

It can be shown that

\begin{equation}\label{reich_5}
A_{1}'= D+ A_2
\end{equation}

Substituting equation (<ref>) in equation (<ref>),

\begin{align}
\bar{z_i}'Az_{i} &=(1-\bar{\mu_i})\bar{z_{i}}'(D+A_2) z_i \nonumber \\
&= (1-\bar{\mu_i})\bar{z_{i}}'Dz_i +
(1-\bar{\mu_i})\bar{z_{i}}'A_2z_i = (1-\bar{\mu_i})\bar{z_{i}}'Dz_i +
\frac{(1-\bar{\mu_i})}{1-\mu_i}\bar{z_i}'Az_i \label{reich_6}
\end{align}

Simplifying equation (<ref>),

\begin{equation}\label{reich_7}
(1-\bar{\mu_i} \mu_i)\bar{z_i}'Az_i = (1-\bar{\mu_i})(1-\mu_i)\bar{z_i}'Dz_i
\end{equation}

But, $\bar{\mu_i}\mu_i= |\mu_i|^2$ and
$(1-\bar{\mu_i})(1-\mu_i)= |1-\mu_i|^2$
Hence, equation (<ref>) simplifies to,

\begin{equation}\label{reich_8}
(1-|\mu_i|^2)\bar{z_i}'Az_i = (|1-\mu_i|^2) \bar{z_i}'Dz_i
\end{equation}

Since, the diagonal elements of $A$ is positive, i.e., $A_{ii}>0$, $\bar{z_i}'Dz_i$ is positive definite, i.e., $\bar{z_i}'Dz_i >0$.

Since $\bar{z_i}'Dz_i >0$ and $\bar{z_i}'Az_i >0$, $|\mu_i|<1$.

§.§ Formalization in Coq

The Theorem <ref> is formalized in Coq as follows:
Theorem Reich_sufficiency: forall (n:nat) (A: 'M[R]_n.+1),
(forall i:'I_n.+1,  A i i > 0) ->
(forall i j:'I_n.+1,   A i j = A j i) ->
is_positive_definite A ->
(let S:= oppmx (mulmx (RtoC_mat (invmx (A1 A))) (RtoC_mat (A2 A))) in
(forall i: 'I_n.+1, C_mod (lambda S i) < 1)).
where positive definiteness of a complex matrix $A$ is defined as:

\begin{equation*}
\forall x \in \mathbb{C}^{n \times 1}, Re\; [x^{*}Ax] >0
\end{equation*}

$x^{*}$ is the complex conjugate transpose of vector $x$ and $Re\; [x^{*}Ax]$ is the real part of the complex scalar $x^{*}Ax$. We defined  in Coq as:
Definition is_positive_definite (n:nat) (A: 'M[R]_n.+1):=
forall (x: 'cV[complex R]_n.+1),  x != 0 ->
Re (mulmx (conjugate_transpose x) (mulmx (RtoC_mat A) x) 0 0) >0.

We compared our definition of a positive definite matrix with a related work from Pierre Roux [31]. While their work define positive definiteness for a real matrix, ours define it for a complex matrix. The definitions however are similar.
The hypothesis
forall i:'I_n.+1,  A i i > 0

states that all terms in the main diagonal of A are positive. The hypothesis
forall i j:'I_n.+1,   A i j = A j i
states that the matrix A is symmetric.

We can then apply the theorem  with  to prove convergence of the Gauss–Seidel iteration method. We formalize the convergence of the Gauss–Seidel iteration method in Coq as
Theorem Gauss_Seidel_converges:
forall (n:nat) (A: 'M[R]_n.+1) (b: 'cV[R]_n.+1),
let x:= (invmx A) *m b in
A ∈unitmx ->
(forall i : 'I_n.+1, 0 < A i i) ->
(forall i j : 'I_n.+1, A i j = A j i) ->
is_positive_definite A ->
(forall x0: 'cV[R]_n.+1,
is_lim_seq (fun m:nat =>
vec_norm (addmx (X_m m.+1 x0 b (A1 A) (A2 A)) (oppmx x))) 0

We next demonstrate the convergence of the Gauss–Seidel iteration on the example (<ref>). We choose $N=1$. Thus, we have a symmetric tri-diagonal coefficient matrix of size $3 \times 3$, which we will denote as $A_{GS}$.
To show that iterative system for the system $A_{GS} x = b$ converges, we need to show that $A_{GS}$ is positive definite. In Coq, we prove that $A_{GS}$ is positive definite as the following lemma statement
Lemma Ah_pd: forall (h:R), (0<h)

Proving that $A_{GS}$ is positive definite by definition for a generic $N$ is tedious and does not add much to our line of argument. Hence, we chose to do it for $A_{GS}$ of size $3 \times 3$. One can perform the exercise for any choice of $N$ and get the same result. The statement of convergence of Gauss–Seidel iteration method for the $3 \times 3$ matrix is stated in Coq as
Theorem Gauss_seidel_Ah_converges:
forall (b: 'cV[R]_3) (h:R),
(0 < h)
let A := (Ah 2
let x:= (invmx A) *m b in
forall x0: 'cV[R]_3,
is_lim_seq (fun m:nat =>
vec_norm (addmx (X_m m.+1 x0 b (A1 A) (A2 A)) (oppmx x))) 0

This closes the proof of the convergence of Gauss–Seidel iteration for the example problem.

§ JACOBI ITERATION

The Jacobi method is another instance of iterative methods to solve the system of linear equations
$ Ax = b \quad \text{or} \quad
\sum_{j=1}^{n} A_{ij}x_j = b_i,\; \forall i= 1 .. n$
by successive approximations. We solve for the values $x_i$ at the $m^{th}$ iteration, by keeping the other elements of the vector $x$ fixed at the value obtained after $m-1$ iterations. This gives us the following update formula [5]:

\begin{equation}\label{jacobi_update}
x_{i,m} = \frac{b_i - \sum_{j \neq i}A_{ij}x_{j, m-1}} {A_{ii}}
\end{equation}

Comparing the system of equations (<ref>) with the iterative system (<ref>), we can define the matrix $A_1$ and $A_2$ as [5]:

\begin{equation}
\begin{cases}
A_{ij} & \text{if $i=j$}\\
0 & \text{otherwise}
\end{cases};
\qquad
\begin{cases}
A_{ij} & \text{if $i \neq j$}\\
0 & \text{otherwise}
\end{cases}
\end{equation}

The iterative matrix $S_J$, would then be:

\begin{equation}\label{jacobi_iterate}
S_J \overset{\Delta}{=} -A_{1}^{-1}A_2 =  -D^{-1}(A-D) = I-D^{-1}A
\end{equation}

where $D$ is the diagonal matrix, $L$ is strictly lower triangular matrix and $U$ is strictly upper triangular matrix of $A$.
Therefore, for the Jacobi method, the Theorem <ref> would be stated as:

Let the Jacobi iterative matrix be defined as in (<ref>) for the iterative system (<ref>). The sequence of iterative solutions $\{x_{m}\}$ converges to the direct solution $x$ if and only if the spectral radius of the iterative matrix (<ref>), $\rho(I-D^{-1}A)$ is less than 1.

We will now formally prove the convergence of the Jacobi iteration for the example problem (<ref>). We choose $N=1$, thereby obtaining a $3 \times 3$ matrix system, like we did for the Gauss–Seidel iteration. However, the Jacobi matrix, $A_J$, that we will consider will be a negation of the coefficient matrix, $A$ in (<ref>). We essentially shift the negation to the right hand side of (<ref>). The reason we do this is that, to define the explicit form of eigenvalues of the Jacobi matrix for the example problem (<ref>), we need to take the square root of the off-diagonal elements of $A_J$. In Coq, the square root definition of a real number has a precondition that the real number must be non-negative. Thus, we require the off-diagonal elements of $A_J$ to be non-negative. This does not change the problem in consideration, but transforms it into a form that is easy to formalize in Coq.

The eigenvalues of $I-D^{-1}A_J$ would then be $1+\frac{h^2}{2} \lambda(A_J)$. For a tri-diagonal matrix with $a$ on the lower diagonal, $b$ on the diagonal and $c$ on the upper diagonal entries, $\lambda_i(A_J)$ is defined as

\begin{equation}\label{eigen}
\lambda_i(A_J)= b + 2 \sqrt{ac}\left[\cos{\left(\frac{m\pi}{N+1}\right)}\right], \quad 0 \leq i <N, \; 1\leq M \leq N
\end{equation}

For the matrix $A_J$, $b=-\frac{2}{h^2}$, $a=c=\frac{1}{h^2}$.

\begin{equation}
\lambda_i(A)=-\frac{2}{h^2}+2\frac{1}{h^2}\left[\cos{\left(\frac{m\pi}{N+1}\right)}\right] = \frac{-2}{h^2}\left[1-\cos{\left(\frac{m\pi}{N+1}\right)}\right]
\end{equation}

\begin{equation}
|\lambda_i(S_J)|= \left| 1+\frac{h^2}{2}\lambda_i(A_J) \right|= \left| \cos{\left(\frac{m\pi}{N+1}\right)}
\right|
\end{equation}

Thus, to prove the convergence of the Jacobi iteration for this example, we need to prove that

\begin{equation}\label{jacobi_inst}
\left|\cos{\left(\frac{m\pi}{N+1}\right)} \right| <1
\end{equation}

§.§ Formalization in Coq

In Coq, we formalized the condition (<ref>) as:
Theorem eig_less_than_1:
forall (n:nat) (i: 'I_n.+1) (h:R),
(0 < h)

where the eigenvalues of the iterative matrix $S_J$ is defined as below:
Definition lambda_J (m N:nat) (h:R) :=
(b N h) + 2* sqrtc(p N h)* RtoC (cos((m.+1

In Coq the natural numbers start from 0, hence we increment $m$ and $N$ by one as compared to the formula (<ref>).

Since the above definition of eigenvalues hold for a tri-diagonal matrix, we formally prove that $S_J = I - D^{-1}A_J$, thereby preserving the tri-diagonal structure of $A_J$. This is stated in Coq as:
Lemma S_mat_simp: forall (n:nat) (h:R), (0<h)
S_mat_J n h =  RtoC_mat (addmx 1

where  is a coercion from a real to a complex matrix as discussed earlier. This lemma holds for $0<h$. This condition is required since $h$ denotes the discretization step size i.e. $h = x_{i+1}-x_{i}$.

Since $\lambda_i (S_J)$ is in general complex, to prove (<ref>), we need to prove that the real part of $\lambda_i(S_J)$ is equal to the the eigenvalue of $A_J$. We prove this formally in Coq as the following lemma:
Lemma lambda_simp: forall (m N:nat) (h:R), (0<h)
Re (lambda_J m N h) = (1 + ((h^2)/2) *Lambda_Ah m N h))
is the $m^{th}$ eigenvalue of the matrix $A_J$. The hypothesis $0< N$ is a constraint on the number of internal points in the domain. The equation (<ref>) dictates that there must be at least one point in the interior.

One of the requirements of defining an iterative system is that the matrix $A_1$ be invertible. We prove this formally in Coq for this example as the following lemma statement:
Lemma A1_invertible: forall (n:nat) (h:R), (0<h)

To prove the above lemma, we have to directly prove that $A_1 A_1^{-1} = I$. This is stated in Coq as the following lemma statement:
Lemma invmx_A1_mul_A1_1:forall (n:nat) (h:R),

Since $A_1$ is a diagonal matrix, we can obtain a direct formulation of the inverse of $A_1$. This is stated in Coq as the following lemma statement:
Lemma invmx_A1_J: forall (n:nat) (h:R), (0<h)
invmx (A1_J n h) = ((-(h^+2)/2) *: 1

The lemma  was the last piece in the puzzle to complete the proof of invertibility of $A_1$. We apply the Theorem <ref> to the $3 \times 3$ example to prove convergence of Jacobi iteration matrix on this example.
Theorem Jacobi_converges: forall (b: 'cV[R]_3) (h:R),
(0 < h)
let A := (Ah 2
let x := (invmx A) *m b in
(forall x0: 'cV[R]_3,
is_lim_seq (fun m:nat =>
vec_norm (addmx (X_m m.+1 x0 b (A1_J 2

Here, we instantiate the eigenvalue  with the  using
Hypothesis Lambda_eq: forall (n:nat) (h:R) (i: 'I_n.+1),
let S_mat :=
(oppmx (invmx (A1_J n h) *m A2_J n h)) in
lambda S_mat i = lambda_J i n h.

We further formally prove that  is an eigenvalue of the Jacobi matrix, $A_J$, using the  predicate in .
Lemma lambda_J_is_an_eigenvalue:
forall (h:R),
(0 < h)
let S_mat := RtoC_mat (oppmx (invmx (A1_J 2 h) *m A2_J 2 h)) in
(forall i: 'I_3,
@eigenvalue (complex_fieldType _) 3

This closes the proof of iterative convergence for Jacobi iteration on the example problem.

§ FORMALIZING PROPERTIES OF COMPLEX NUMBERS, COMPLEX MATRICES AND VECTORS

While there exist a rich formalization of real analysis, linear algebra, finite sets and sequences in Coq, generic properties about complex vectors and matrices were lacking. In this section, we will discuss about how we addressed these gaps in our formalization.

§.§ Formalizing properties of complex matrices and vectors

The  theory in the  [28] library in  define complex numbers and basic operations on them. They define complex numbers as a real closed field, thereby allowing us to instantiate a generic field with a complex field. This was useful when we used the  definition from  matrix algebra library and the canonical forms library by Cano et al [13]. However, since the
basic properties like modulus of a complex number, conjugates, properties of complex matrices and vectors were lacking, we added them in our formalization.

We provided the definition of the modulus of a complex number in Section <ref>. We proved some basic properties of the modulus, which we enumerate in Table <ref>.

Formalization of properties of complex modulus in Coq

Mathematical properties     Formalization in Coq

$|| 0 || = 0$        [keywordstyle=,language=Coq]
Lemma C_mod_0: C_mod 0 = 0

$0 \leq ||x||$
Lemma C_mod_ge_0:
forall (x: complex R), (0<= C_mod x)

$||x y|| = ||x||\; ||y||$
Lemma C_mod_prod: forall (x y: complex R),
C_mod (x * y) = C_mod x * C_mod y.

$||\frac{x}{y}|| =  \frac{||x||}{||y||}, y \neq 0$
Lemma C_mod_div: forall (x y: complex R),
y <> 0 ->
C_mod (x / y) = (C_mod x) / (C_mod y).

$||x|| \neq 0,\quad \text{if. } x \neq 0$
Lemma C_mod_not_zero: forall (x: complex R),
x <> 0 -> C_mod x <> 0.

$||1|| = 1$
Lemma C_mod_1: C_mod 1 = 1.

$||x^n|| = ||x||^n$
Lemma C_mod_pow: forall (x: complex R) (n:nat),
C_mod (x^+ n) = (C_mod x)^+n.

$||x + y|| \leq ||x|| + ||y||$
Lemma C_mod_add_leq : forall (a b: complex R),
C_mod (a + b) <= C_mod a + C_mod b.

$||\frac{1}{x}|| = \frac{1}{||x||}, \quad \text{if  } x \neq 0$
Lemma C_mod_inv : forall x : complex R,
x <> 0 -> C_mod (invc x) = Rinv (C_mod x).

$||xy||^2 = ||x||^2 ||y||^2$
Lemma C_mod_sqr: forall (x y : complex R),
Rsqr (C_mod (x * y)) =
(Rsqr (C_mod x)) * (Rsqr (C_mod y)).

$||-x|| = ||x||$
Lemma C_mod_minus_x: forall (x: complex R),
C_mod (-x) = C_mod x.

$||\sum_{j=0}^n u(j)|| \leq
\sum_{j=0}^n ||u(j)||
Lemma C_mod_sum_rel:
forall (n:nat) (u : 'I_n.+1 -> (complex R)),
(C_mod ([+

The  theory in  defines the conjugate $\bar{x}$ of a complex number $x$ as
Definition conjc R : ringType (x : R[i]) := let: a +i* b := x in a -i* b.

The Table <ref> lists the missing formalization of complex conjugates that we added for this formalization.

Formalization of properties of complex conjugates in Coq

Mathematical properties     Formalization in Coq

$\overline{xy} = \bar{x} \; \bar{y}$        [keywordstyle=,language=Coq]
Lemma Cconj_prod: forall (x y: complex R),
conjc (x*y)

$\overline{x + y} = \bar{x} + \bar{y} $
Lemma Cconj_add: forall (x y: complex R),
conjc (x+y) = conjc x + conjc y.

$||x|| = ||\bar{x}||$
Lemma Cconjc_mod: forall (a: complex R),
C_mod a = C_mod (conjc a).

$x = \bar{\bar{x}}$
Lemma conj_of_conj_C: forall (x: complex R),
x = conjc (conjc x).

$\bar{x} x = ||x||^2$
Lemma conj_prod: forall (x:complex R),
((conjc x)*x)

$Re[x] + Re[\bar{x}] = 2 Re[x]$
Lemma Re_conjc_add: forall (x: complex R),
Re x + Re (conjc x) = 2 * (Re x).

$\overline{\sum_{j=0}^n f(i)} = \sum_{j=0}^n \overline{f(i)}$
Lemma Cconj_sum: forall (p:nat) (x: 'I_p -> complex R),
conjc ([+

Here,  is a coercion from reals to complex.

We define the conjugate transpose of a complex matrix in Coq as
Definition conjugate_transpose (m n:nat) (A: 'M[complex R]_(m,n)):=
transpose_C (conjugate A).

where  is the transpose of a complex matrix, which we define in Coq as
Definition transpose_C (m n:nat) (A: 'M[complex R]_(m,n)):=
_(i<n,j<m) A j i.

and  is the conjugate of a rectangular matrix. In Coq,
Definition conjugate (m n:nat) (A: 'M[complex R]_(m,n)):=
_(i<m,j<n) conjc (A i j).

The lemma
Lemma conj_scal_mat_mul:
forall (m n : nat) (l:complex R) (x: 'M[complex R]_(m,n)),
conjugate_transpose (scal_mat_C l x) =
scal_mat_C (conjc l) (conjugate_transpose x).

proves the scaling property of a complex matrix, $A$

\begin{equation*}
(\overline{lA})^T = \bar{l}(\bar{A})^T
\end{equation*}

Lemma conj_matrix_mul :
forall (m n p:nat) (A: 'M[complex R]_(m,p)) (B: 'M[complex R]_(p,n)),
conjugate_transpose (mulmx A B) = mulmx
(conjugate_transpose B) (conjugate_transpose A).

The lemma  states that the conjugate transpose of the product of matrices $A$ and $B$ equals the product of conjugate transpose of the matrices, i.e., $\overline{AB} = \bar{B}\bar{A}$.
Lemma conj_of_conj: forall (m n:nat) (x: 'M[complex R]_(m,n)),
conjugate_transpose (conjugate_transpose x) = x.

The lemma  states the idempotent property of the conjugate transpose of a complex vector, $v$, i.e., $\overline{\overline{v}} = v$.

§.§ Formalization of vector and matrix norms

Another missing piece in the existing formalization was norm of a vector and a matrix.
We provided a formal definition of the norm of a complex matrix and vector in Section <ref>.

In Table <ref> and Table <ref>, we enumerate the properties of matrix and vector norms that we formalized.

Formalization of properties of matrix norm in Coq

Mathematical properties     Formalization in Coq

$0 \leq ||A||_i$        [keywordstyle=,language=Coq]
Lemma matrix_norm_ge_0:
forall (n:nat) (A: 'M[complex R]_n.+1),
(0 <= matrix_norm A)

$||Ax|| \leq ||A||_i ||x||, \quad x \neq 0 $[1]
Lemma matrix_norm_compat:
forall (n:nat) (x: 'cV[complex R]_n.+1)
(A: 'M[complex R]_n.+1),
x != 0 ->
vec_norm_C (mulmx A x) <=
((matrix_norm A) * vec_norm_C x)

$||AB||_i \leq ||A||_i ||B||_i$
Lemma matrix_norm_prod:
forall (n:nat) (A B: 'M[complex R]_n.+1),
(matrix_norm (A *m B) <=
(matrix_norm A) * (matrix_norm B))

$0 \leq ||A||_i \leq ||A||_F$[2]
Lemma mat_2_norm_F_norm_compat:
forall (n:nat) (A: 'M[complex R]_n.+1),
(0 <= matrix_norm A <= mat_norm A)

[1]Here, $x$ is a vector and the relation proves compatibilty relation between a matrix norm and its induced vector norm.
[2]Here, $||A||_F$ is the Frobenius norm and we prove that the 2-norm of a matrix is bounded above by the Frobenius matrix norm.
The Frobenius norm of a matrix is defined as

\begin{equation*}
||A||_F = \sqrt{\sum_{j=1}^n \sum_{j=1}^n |A_{ij}|^2}
\end{equation*}

An important point to note here is that since we are using the  definition of an extended real line, , coercion of a quantity of type  to real requires us to prove finiteness of that quantity. We therefore have to prove that the matrix norm is finite, which we state as the following lemma in Coq
Lemma matrix_norm_is_finite: forall (n:nat) (A: 'M[complex R]_n.+1),
is_finite (matrix_norm A).

Formalization of properties of vector norm in Coq

Mathematical properties     Formalization in Coq

$0 \leq ||v||$        [keywordstyle=,language=Coq]
Lemma vec_norm_C_ge_0:
forall (n:nat) (v: 'cV[complex R]_n.+1),
(0<= vec_norm_C v)

$||av|| = |a|\; ||v||, \quad \text{$a$ is scalar} $
Lemma ei_vec_ei_compat: forall (n:nat)
(x:complex R) (v: 'cV[complex R]_n.+1),
vec_norm_C (scal_vec_C x v) =
C_mod x * vec_norm_C v.

$||v_1 + v_2|| \leq ||v_1|| + ||v_2||$
Lemma vec_norm_add_le:
forall (n:nat) (v1 v2 : 'cV[complex R]_n.+1),
vec_norm_C (v1 + v2) <=
vec_norm_C v1 + vec_norm_C v2.

$v \neq 0 \implies ||v|| \neq 0$[1]
Lemma non_zero_vec_norm: forall (n:nat)
(v: 'cV[complex R]_n.+1),
vec_not_zero v -> (vec_norm_C v <> 0)

[1] is Coq's definition of $v \neq 0$.

In Coq, we define the Frobenius matrix norm as
Definition mat_norm (n:nat) (A: 'M[complex R]_n.+1) : R:=
sqrt (∑_i (∑_j (Rsqr (C_mod (A i j))))).

§ CONCLUSION AND FUTURE WORK

In this work we formalized the iterative convergence error.
We formalized the necessary and sufficient conditions for the convergence of a sequence of iterative solutions $\{x_m\}$ to the direct solution $x$ of the original linear system $Ax=b$. To prove convergence, one needs to construct an iterative matrix $S\overset{\Delta}{=}A_{1}^{-1}A_2$ and prove that the magnitude of its eigenvalues is less than 1. The definitions of matrices $A_1$ and $A_2$ depend on the choice of iterative methods. In this paper, we provide a detailed proof, especially for the convergence of the Jordan block diagonal matrix for the magnitude of eigenvalues less than 1. In our formal proof we clarify some details regarding this convergence.
In this work, we particularly looked at the Gauss–Seidel iterative process. We then formalized the Reich theorem which guarantees iterative convergence if the original matrix $A$ is positive definite. We then instantiated the Gauss–Seidel iteration with an example and applied the iterative convergence and the Reich theorem to prove convergence on this example. We also proved the convergence of the Jacobi iterative process for the same example problem. During our formalization, we develop a library in Coq to deal with complex vectors and matrices. We defined absolute values of complex numbers, common properties of complex conjugates and operations on conjugate matrices and vectors. This development leverages the existing formalization [28, 29] of complex numbers and matrices in . The overall length of the Coq code and proofs is about 8.5k lines of code. It took us about 8 person-months of full time work for the entire formalization.

As an extension to this work, we are working on incorporating the effect of floating-point arithmetic on the iterative convergence. We plan on using the tools VCFloat [32] and VST [33] to generate proof obligations in Coq to prove that the actual implementation of iterative algorithms obey the high level specifications.
Another possible direction would be to extend the formalization to deal with advanced iterative process like Krylov space methods and other complicated differential equation solvers.
In this work we looked at the linear differential equations. In future work, we would like to also look at non-linear conditions and systems. Most physical systems behave non-linearly. Therefore, it would be interesting to go through the process of linearizing them about the optimal points in a formal setting and formalize conditions for existence and uniqueness of their solutions.

One of the overarching goals that we would like to achieve is to verify the algorithms used by widely used ordinary differential equation solvers like SUNDIALS [7], ODEPACK [6] etc. Most of the high performance computing and cyber-physical applications rely on these common purpose solvers for computing solutions of the differential equations. Therefore, verifying their algorithms will provide rigorous formal guarantees and expose possible bugs in the implementation of the algorithms.

§.§.§ Acknowledgments

We thank Yves Bertot for porting the canonical forms library to newer version of Coq, and helping us with understanding the canonical forms and the mathcomp library.

#1ISBN #1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1<https://doi.org/#1>et al.#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1<><#>1#1#1#1#1#1#1#1#1#1#1#1#1#1PreBibitemsHook

[1]

Batchelor, G.K.:
An Introduction to Fluid Dynamics,
2nd pbk. ed. edn.
Cambridge mathematical library.
Cambridge University Press,
Cambridge, U.K., New York, NY

[2]

Reid, W.T.:
Riccati Differential Equations.
Mathematics in science and engineering, v. 86.
Academic Press,
New York

[3]

Wilmott, P.,
Howison, S.,
Dewynne, J.:
The Mathematics of Financial Derivatives: A Student Introduction.
Cambridge University Press,

[4]

The Prize in Economic Sciences 1997 - Press release.
(Accessed on 09/22/2021)

[5]

Saad, Y.:
Iterative Methods for Sparse Linear Systems,
2nd ed. edn.
Society for Industrial and Applied Mathematics (SIAM),
Philadelphia, Pa.

[6]

Hindmarsh, A.C.:
Odepack, a systematized collection of ode solvers.
Scientific computing,

[7]

Hindmarsh, A.C.,
Brown, P.N.,
Grant, K.E.,
Lee, S.L.,
Serban, R.,
Shumaker, D.E.,
Woodward, C.S.:
SUNDIALS: Suite of nonlinear and differential/algebraic equation
ACM Transactions on Mathematical Software (TOMS)

[8]

O’Connor, R.:
Certified exact transcendental real number computation in coq.
In: International Conference on Theorem Proving in Higher Order
pp. 246–261

[9]

Boldo, S.,
Lelay, C.,
Melquiond, G.:
Coquelicot: A user-friendly library of real analysis for Coq.
Mathematics in Computer Science

[10]

Garillot, F.,
Gonthier, G.,
Mahboubi, A.,
Rideau, L.:
Packaging mathematical structures.
In: International Conference on Theorem Proving in Higher Order
pp. 327–342

[11]

Martin-Dorel, É.,
Rideau, L.,
Théry, L.,
Mayero, M.,
Pasca, I.:
Certified, efficient and sharp univariate taylor models in coq.
In: 2013 15th International Symposium on Symbolic and Numeric
Algorithms for Scientific Computing,
pp. 193–200

[12]

Pasca, I.:
Formal verification for numerical methods.
PhD thesis,
Université Nice Sophia Antipolis

[13]

Cano, G.,
Dénès, M.:
Matrices à blocs et en forme canonique.
In: Pous, D.,
Tasson, C. (eds.)
JFLA - Journées Francophones des Langages Applicatifs,
Aussois, France
Damien Pous and Christine Tasson.

[14]

Thiemann, R.:
A perron–frobenius theorem for deciding matrix growth.
Journal of Logical and Algebraic Methods in Programming,

[15]

Bréhard, F.:
Numerical computation certified in functional spaces: A trilogue between
rigorous polynomial approximations, symbolic computation and formal proof.
PhD thesis

[16]

Boldo, S.,
Clément, F.,
Filliâtre, J.-C.,
Mayero, M.,
Melquiond, G.,
Weis, P.:
Formal proof of a wave equation resolution scheme: the method error.
In: International Conference on Interactive Theorem Proving,
pp. 147–162

[17]

Boldo, S.,
Clément, F.,
Filliâtre, J.-C.,
Mayero, M.,
Melquiond, G.,
Weis, P.:
Trusting computations: a mechanized proof from partial differential
equations to actual program.
Computers & Mathematics with Applications

[18]

Boldo, S.,
Clément, F.,
Filliâtre, J.-C.,
Mayero, M.,
Melquiond, G.,
Weis, P.:
Wave equation numerical resolution: a comprehensive mechanized proof
of a c program.
Journal of Automated Reasoning

[19]

Tekriwal, M.,
Duraisamy, K.,
Jeannin, J.-B.:
A formal proof of the lax equivalence theorem for finite difference
In: Dutle, A.,
Moscato, M.M.,
Titolo, L.,
Muñoz, C.A.,
Perez, I. (eds.)
NASA Formal Methods,
pp. 322–339.

[20]

Immler, F.,
Hölzl, J.:
Numerical analysis of ordinary differential equations in
In: International Conference on Interactive Theorem Proving,
pp. 377–392

[21]

Immler, F.:
A verified ode solver and smale's 14th problem.
Technische Universität München,

[22]

Immler, F.,
Traut, C.:
The flow of odes.
In: International Conference on Interactive Theorem Proving,
pp. 184–199

[23]

Immler, F.,
Traut, C.:
The flow of odes: Formalization of variational equation and
poincaré map.
Journal of Automated Reasoning

[24]

Immler, F.:
Formally verified computation of enclosures of solutions of ordinary
differential equations.
In: Badger, J.M.,
Rozier, K.Y. (eds.)
NASA Formal Methods,
pp. 113–127.

[25]

Mahboubi, A.,
Tassi, E.:
Mathematical components

[26]

Sohrab, H.H.:
Basic Real Analysis.

[27]

The Ratio Test for Sequence Convergence - Mathonline.
(Accessed on 09/21/2021)

[28]

Cohen, C.:
Construction of real algebraic numbers in Coq.
In: Beringer, L.,
Felty, A. (eds.)
ITP - 3rd International Conference on Interactive Theorem Proving -
Princeton, United States

[29]

Cohen, C.,
Mahboubi, A.:
Formal proofs in real algebraic geometry: from ordered fields to
quantifier elimination.
Logical Methods in Computer Science

[30]

Reich, E.:
On the convergence of the classical iterative method of solving linear
simultaneous equations.
The Annals of Mathematical Statistics

[31]

Roux, P.:
Formal proofs of rounding error bounds.
Journal of Automated Reasoning

[32]

Ramananandro, T.,
Mountcastle, P.,
Meister, B.,
Lethin, R.:
A unified coq framework for verifying c programs with floating-point
In: Proceedings of the 5th ACM SIGPLAN Conference on Certified
Programs and Proofs.
CPP 2016,
pp. 15–26.
Association for Computing Machinery,
New York, NY, USA

[33]

Appel, A.W.:
Verified software toolchain.
In: Barthe, G. (ed.)
Programming Languages and Systems,
pp. 1–17.
Berlin, Heidelberg

§ INFORMAL PROOF OF THE RATIO TEST

Here, we provide an informal proof [27] of the Theorem (<ref>).

Let $(a_n)$ be a sequence of real numbers such that $\lim_{n \to \infty} \frac{a_{n+1}}{a_n} = L$. Since $(a_n)$ is a sequence of positive numbers, $0 \leq L$. By the density of Real Numbers theorem, there exists a real $r$, such that $L < r < 1$.

Unfolding the definition of $\lim_{n \to \infty} \frac{a_{n+1}}{a_n} = L$, for $\epsilon = r- L > 0 $, $\exists N \in \mathbb{N}$, such that forall $n \geq N$, $\left| \frac{a_{n+1}}{a_n} - L \right| < \epsilon = r - L$,  which implies that,

\begin{equation*}
-\epsilon < \frac{a_{n+1}}{a_n} - L < \epsilon
\implies L- \epsilon < \frac{a_{n+1}}{a_n} < L + \epsilon
\implies \frac{a_{n+1}}{a_n} < L + (r-L) = r
\end{equation*}

Therefore, $\forall n, \; n \geq N$, $a_{n+1} < a_n r$.

This implies, $a_n r < a_{n-1}r^2$, and $a_{n-2}r^2 < a_{n-3}r^3$, ...,
$a_{N+1}r^{n-N} < a_N r^{n-N+1}$

We thus obtain the following inequality:

\begin{equation*}
a_{n+1} < a_n r < a_{n-1}r^2 < \hdots < a_{N+1}r^{n-N} < a_N r^{n-N+1}
\end{equation*}

Let $K= \frac{a_N}{r^N}$. Therefore $a_N r^{n-N+1} = K r^{n+1}$. Thus for $n \geq N$, we have that $a_n r < K r^{n+1}$, or rather $a_n < K r^n$. Since $0 < r < 1$, we have that $\lim_{n \to \infty} r^n = 0$. Then by the squeeze lemma, $\lim_{n \to \infty} a_n = 0.$

§ INFORMAL PROOF OF (<REF>)

To prove that the inequality (<ref>) holds, unfold the definition of ${m \choose k}$. Therefore,

\begin{align}
{m \choose k} &= \frac{m!}{(m-k)! k!} \nonumber \\
&= \frac{(m-k+1)(m-k+2)...(m-1)m}{k!}\label{choose_num}
\end{align}

Since $(m-k+1) \leq m$,  $(m-k+2) \leq m$, and so on,
the product of terms in the numerator of (<ref>) is bounded by $m^k$. Hence,

\begin{equation*}
{m \choose k} \leq
\frac{m^k}{k!}
\end{equation*}

§ RESULTS

Sample body text. Sample body text. Sample body text. Sample body text. Sample body text. Sample body text. Sample body text. Sample body text.

§ THIS IS AN EXAMPLE FOR FIRST LEVEL HEAD—SECTION HEAD

§.§ This is an example for second level head—subsection head

§.§.§ This is an example for third level head—subsubsection head

Sample body text. Sample body text. Sample body text. Sample body text. Sample body text. Sample body text. Sample body text. Sample body text.

§ EQUATIONS

Equations in  can either be inline or on-a-line by itself (“display equations”). For
inline equations use the  commands. E.g.: The equation
$H\psi = E \psi$ is written via the command .

For display equations (with auto generated equation numbers)
one can use the equation or align environments:

\begin{equation}
\|\tilde{X}(k)\|^2 \leq\frac{\sum\limits_{i=1}^{p}\left\|\tilde{Y}_i(k)\right\|^2+\sum\limits_{j=1}^{q}\left\|\tilde{Z}_j(k)\right\|^2 }{p+q}.\label{eq1}
\end{equation}

\begin{align}
D_\mu &=  \partial_\mu - ig \frac{\lambda^a}{2} A^a_\mu \nonumber \\
F^a_{\mu\nu} &= \partial_\mu A^a_\nu - \partial_\nu A^a_\mu + g f^{abc} A^b_\mu A^a_\nu \label{eq2}
\end{align}

Notice the use of  in the align environment at the end
of each line, except the last, so as not to produce equation numbers on
lines where no equation numbers are required. The  command
should only be used at the last line of an align environment where
is not used.

\begin{equation}
Y_\infty = \left( \frac{m}{\textrm{GeV}} \right)^{-3}
\left[ 1 + \frac{3 \ln(m/\textrm{GeV})}{15}
+ \frac{\ln(c_2/5)}{15} \right]
\end{equation}

The class file also supports the use of ,  and
commands. As such ,
and  produces $\mathbb{R}$, $\mathscr{R}$ and $\mathcal{R}$
respectively (refer Subsubsection <ref>).

§ TABLES

Tables can be inserted via the normal table and tabular environment. To put
footnotes inside tables you should use  tag.
The footnote appears just below the table itself (refer Tables <ref> and <ref>).
For the corresponding footnotemark use

Caption text

Column 1     Column 2      Column 3     Column 4

row 1        data 1       data 2      data 3

row 2        data 4       data 5[1]      data 6

row 3        data 7       data 8      data 9[2]

Source: This is an example of table footnote. This is an example of table footnote.
[1]Example for a first table footnote. This is an example of table footnote.
[2]Example for a second table footnote. This is an example of table footnote.

The input format for the above table is as follows:

Example of a lengthy table which is set to full textwidth

3@c@Element 1[1]     3@c@Element 2[2]
Project     Energy     $\sigma_{calc}$     $\sigma_{expt}$     Energy     $\sigma_{calc}$     $\sigma_{expt}$

Element 3      990 A     1168     $1547\pm12$     780 A     1166     $1239\pm100$

Element 4      500 A     961      $922\pm10$      900 A     1268     $1092\pm40$

Note: This is an example of table footnote. This is an example of table footnote this is an example of table footnote this is an example of table footnote this is an example of table footnote.
[1]Example for a first table footnote.
[2]Example for a second table footnote.

In case of double column layout, tables which do not fit in single column width should be set to full text width. For this, you need to use    instead of    environment. Lengthy tables which do not fit in textwidth should be set as rotated table. For this, you need to use    instead of    environment. This environment puts tables rotated to single column width. For tables rotated to double column width, use   .

Tables which are too long to fit, should be written using the “sidewaystable” environment as shown here

3@c@Element 1[1]    3@c@Element[2]
Projectile     Energy	    $\sigma_{calc}$     $\sigma_{expt}$     Energy     $\sigma_{calc}$     $\sigma_{expt}$

Element 3     990 A     1168     $1547\pm12$     780 A     1166     $1239\pm100$

Element 4     500 A     961      $922\pm10$      900 A     1268     $1092\pm40$

Element 5     990 A     1168     $1547\pm12$     780 A     1166     $1239\pm100$

Element 6     500 A     961      $922\pm10$      900 A     1268     $1092\pm40$

Note: This is an example of table footnote this is an example of table footnote this is an example of table footnote this is an example of table footnote this is an example of table footnote.
[1]This is an example of table footnote.

§ FIGURES

As per the  standards you need to use eps images for  compilation and  images for  compilation. This is one of the major difference between  and . Each image should be from a single input .eps/vector image file. Avoid using subfigures. The command for inserting images for  and  can be generalized. The package used to insert images in  is the graphicx package. Figures can be inserted via the normal figure environment as shown in the below example:

This is a widefig. This is an example of long caption this is an example of long caption  this is an example of long caption this is an example of long caption

In case of double column layout, the above format puts figure captions/images to single column width. To get spanned images, we need to provide   .

For sample purpose, we have included the width of images in the optional argument of  tag. Please ignore this.

§ ALGORITHMS, PROGRAM CODES AND LISTINGS

Packages ,  and  are used for setting algorithms in  using the format:

You may refer above listed package documentations for more details before setting  environment. For program codes, the “program” package is required and the command to be used is   . A fast exponentiation procedure:

i:=1 10 1 |expt|(2,i);
Comments will be set flush to the right margin
|expt|(x,n) z:=1;
n=0 ;
|odd|(n) ;
This is a comment statement;
n:=n/2; x:=x*x ;
{ n>0 };
n:=n-1; z:=z*x ;

Calculate $y = x^n$
[1]
$n \geq 0 \vee x \neq 0$
$y = x^n$
$y \Leftarrow 1$
$n < 0$
$X \Leftarrow 1 / x$
$N \Leftarrow -n$
$X \Leftarrow x$
$N \Leftarrow n$
$N \neq 0$
$N$ is even
$X \Leftarrow X \times X$
$N \Leftarrow N / 2$
[$N$ is odd]
$y \Leftarrow y \times X$
$N \Leftarrow N - 1$

Similarly, for , use the  package.    is used to set environments similar to  environment. Refer to the  package documentation for more details.

for i:=maxint to 0 do
do nothing
Write('Case insensitive ');
Write('Pascal keywords.');

§ CROSS REFERENCING

Environments such as figure, table, equation and align can have a label
declared via the  command. For figures and table
environments use the  command inside or just
below the  command. You can then use the
command to cross-reference them. As an example, consider
the label declared for Figure <ref> which is
. To cross-reference it, use the command
, for which it comes up as
“Figure <ref>”.

To reference line numbers in an algorithm, consider the label declared for the line number 2 of Algorithm <ref> is . To cross-reference it, use the command  for which it comes up as line <ref> of Algorithm <ref>.

§.§ Details on reference citations

Standard  permits only numerical citations. To support both numerical and author-year citations this template uses   package. For style guidance please refer to the template user manual.

Here is an example for : <cit.>. Another example for : <cit.>. For author-year citation mode,  prints Jones et al. (1990) and  prints (Jones et al., 1990).

All cited bib entries are printed at the end of this article: <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.> and <cit.>.

§ EXAMPLES FOR THEOREM LIKE ENVIRONMENTS

For theorem like environments, we require  package. There are three types of predefined theorem styles exists—,  and

Numbered, theorem head in bold font and theorem text in italic style
Numbered, theorem head in roman font and theorem text in italic style
Numbered, theorem head in bold font and theorem text in roman style

For mathematics journals, theorem styles can be included as shown in the following examples:

Example theorem text. Example theorem text. Example theorem text. Example theorem text. Example theorem text.
Example theorem text. Example theorem text. Example theorem text. Example theorem text. Example theorem text.
Example theorem text.

Sample body text. Sample body text. Sample body text. Sample body text. Sample body text. Sample body text. Sample body text. Sample body text.

Example proposition text. Example proposition text. Example proposition text. Example proposition text. Example proposition text.
Example proposition text. Example proposition text. Example proposition text. Example proposition text. Example proposition text.

Sample body text. Sample body text. Sample body text. Sample body text. Sample body text. Sample body text. Sample body text. Sample body text.

Phasellus adipiscing semper elit. Proin fermentum massa
ac quam. Sed diam turpis, molestie vitae, placerat a, molestie nec, leo. Maecenas lacinia. Nam ipsum ligula, eleifend
at, accumsan nec, suscipit a, ipsum. Morbi blandit ligula feugiat magna. Nunc eleifend consequat lorem.

Sample body text. Sample body text. Sample body text. Sample body text. Sample body text. Sample body text. Sample body text. Sample body text.

Phasellus adipiscing semper elit. Proin fermentum massa
ac quam. Sed diam turpis, molestie vitae, placerat a, molestie nec, leo. Maecenas lacinia. Nam ipsum ligula, eleifend
at, accumsan nec, suscipit a, ipsum. Morbi blandit ligula feugiat magna. Nunc eleifend consequat lorem.

Sample body text. Sample body text. Sample body text. Sample body text. Sample body text. Sample body text. Sample body text. Sample body text.

Example definition text. Example definition text. Example definition text. Example definition text. Example definition text. Example definition text. Example definition text. Example definition text.

Additionally a predefined “proof” environment is available:   . This prints a “Proof” head in italic font style and the “body text” in roman font style with an open square at the end of each proof environment.

Example for proof text. Example for proof text. Example for proof text. Example for proof text. Example for proof text. Example for proof text. Example for proof text. Example for proof text. Example for proof text. Example for proof text.

Sample body text. Sample body text. Sample body text. Sample body text. Sample body text. Sample body text. Sample body text. Sample body text.

Example for proof text. Example for proof text. Example for proof text. Example for proof text. Example for proof text. Example for proof text. Example for proof text. Example for proof text. Example for proof text. Example for proof text.

For a quote environment, use

Quoted text example. Aliquam porttitor quam a lacus. Praesent vel arcu ut tortor cursus volutpat. In vitae pede quis diam bibendum placerat. Fusce elementum
convallis neque. Sed dolor orci, scelerisque ac, dapibus nec, ultricies ut, mi. Duis nec dui quis leo sagittis commodo.

Sample body text. Sample body text. Sample body text. Sample body text. Sample body text (refer Figure <ref>). Sample body text. Sample body text. Sample body text (refer Table <ref>).

§ METHODS

Topical subheadings are allowed. Authors must ensure that their Methods section includes adequate experimental and characterization data necessary for others in the field to reproduce their work. Authors are encouraged to include RIIDs where appropriate.

Ethical approval declarations (only required where applicable) Any article reporting experiment/s carried out on (i) live vertebrate (or higher invertebrates), (ii) humans or (iii) human samples must include an unambiguous statement within the methods section that meets the following requirements:

* Approval: a statement which confirms that all experimental protocols were approved by a named institutional and/or licensing committee. Please identify the approving body in the methods section

* Accordance: a statement explicitly saying that the methods were carried out in accordance with the relevant guidelines and regulations

* Informed consent (for experiments involving humans or human tissue samples): include a statement confirming that informed consent was obtained from all participants and/or their legal guardian/s

If your manuscript includes potentially identifying patient/participant information, or if it describes human transplantation research, or if it reports results of a clinical trial then  additional information will be required. Please visit (<https://www.nature.com/nature-research/editorial-policies>) for Nature Portfolio journals, (<https://www.springer.com/gp/authors-editors/journal-author/journal-author-helpdesk/publishing-ethics/14214>) for Springer Nature journals, or (<https://www.biomedcentral.com/getpublished/editorial-policies#ethics+and+consent>) for BMC.

§ DISCUSSION

Discussions should be brief and focused. In some disciplines use of Discussion or `Conclusion' is interchangeable. It is not mandatory to use both. Some journals prefer a section `Results and Discussion' followed by a section `Conclusion'. Please refer to Journal-level guidance for any specific requirements.

§ CONCLUSION

Conclusions may be used to restate your hypothesis or research question, restate your major findings, explain the relevance and the added value of your work, highlight any limitations of your study, describe future directions for research and recommendations.

In some disciplines use of Discussion or 'Conclusion' is interchangeable. It is not mandatory to use both. Please refer to Journal-level guidance for any specific requirements.

Supplementary information

If your article has accompanying supplementary file/s please state so here.

Authors reporting data from electrophoretic gels and blots should supply the full unprocessed scans for key as part of their Supplementary information. This may be requested by the editorial team/s if it is missing.

Please refer to Journal-level guidance for any specific requirements.

Acknowledgments are not compulsory. Where included they should be brief. Grant or contribution numbers may be acknowledged.

Please refer to Journal-level guidance for any specific requirements.

§ DECLARATIONS

Some journals require declarations to be submitted in a standardised format. Please check the Instructions for Authors of the journal to which you are submitting to see if you need to complete this section. If yes, your manuscript must contain the following sections under the heading `Declarations':

* Funding

* Conflict of interest/Competing interests (check journal-specific guidelines for which heading to use)

* Ethics approval

* Consent to participate

* Consent for publication

* Availability of data and materials

* Code availability

* Authors' contributions

If any of the sections are not relevant to your manuscript, please include the heading and write `Not applicable' for that section.

Editorial Policies for:

Springer journals and proceedings: <https://www.springer.com/gp/editorial-policies>

Nature Portfolio journals: <https://www.nature.com/nature-research/editorial-policies>

Scientific Reports: <https://www.nature.com/srep/journal-policies/editorial-policies>

BMC journals: <https://www.biomedcentral.com/getpublished/editorial-policies>

§ SECTION TITLE OF FIRST APPENDIX

An appendix contains supplementary information that is not an essential part of the text itself but which may be helpful in providing a more comprehensive understanding of the research problem or it is information that is too cumbersome to be included in the body of the paper.

#1ISBN #1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1<https://doi.org/#1>et al.#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1<><#>1#1#1#1#1#1#1#1#1#1#1#1#1#1PreBibitemsHook

[1]

Batchelor, G.K.:
An Introduction to Fluid Dynamics,
2nd pbk. ed. edn.
Cambridge mathematical library.
Cambridge University Press,
Cambridge, U.K., New York, NY

[2]

Reid, W.T.:
Riccati Differential Equations.
Mathematics in science and engineering, v. 86.
Academic Press,
New York

[3]

Wilmott, P.,
Howison, S.,
Dewynne, J.:
The Mathematics of Financial Derivatives: A Student Introduction.
Cambridge University Press,

[4]

The Prize in Economic Sciences 1997 - Press release.
(Accessed on 09/22/2021)

[5]

Saad, Y.:
Iterative Methods for Sparse Linear Systems,
2nd ed. edn.
Society for Industrial and Applied Mathematics (SIAM),
Philadelphia, Pa.

[6]

Hindmarsh, A.C.:
Odepack, a systematized collection of ode solvers.
Scientific computing,

[7]

Hindmarsh, A.C.,
Brown, P.N.,
Grant, K.E.,
Lee, S.L.,
Serban, R.,
Shumaker, D.E.,
Woodward, C.S.:
SUNDIALS: Suite of nonlinear and differential/algebraic equation
ACM Transactions on Mathematical Software (TOMS)

[8]

O’Connor, R.:
Certified exact transcendental real number computation in coq.
In: International Conference on Theorem Proving in Higher Order
pp. 246–261

[9]

Boldo, S.,
Lelay, C.,
Melquiond, G.:
Coquelicot: A user-friendly library of real analysis for Coq.
Mathematics in Computer Science

[10]

Garillot, F.,
Gonthier, G.,
Mahboubi, A.,
Rideau, L.:
Packaging mathematical structures.
In: International Conference on Theorem Proving in Higher Order
pp. 327–342

[11]

Martin-Dorel, É.,
Rideau, L.,
Théry, L.,
Mayero, M.,
Pasca, I.:
Certified, efficient and sharp univariate taylor models in coq.
In: 2013 15th International Symposium on Symbolic and Numeric
Algorithms for Scientific Computing,
pp. 193–200

[12]

Pasca, I.:
Formal verification for numerical methods.
PhD thesis,
Université Nice Sophia Antipolis

[13]

Cano, G.,
Dénès, M.:
Matrices à blocs et en forme canonique.
In: Pous, D.,
Tasson, C. (eds.)
JFLA - Journées Francophones des Langages Applicatifs,
Aussois, France
Damien Pous and Christine Tasson.

[14]

Thiemann, R.:
A perron–frobenius theorem for deciding matrix growth.
Journal of Logical and Algebraic Methods in Programming,

[15]

Bréhard, F.:
Numerical computation certified in functional spaces: A trilogue between
rigorous polynomial approximations, symbolic computation and formal proof.
PhD thesis

[16]

Boldo, S.,
Clément, F.,
Filliâtre, J.-C.,
Mayero, M.,
Melquiond, G.,
Weis, P.:
Formal proof of a wave equation resolution scheme: the method error.
In: International Conference on Interactive Theorem Proving,
pp. 147–162

[17]

Boldo, S.,
Clément, F.,
Filliâtre, J.-C.,
Mayero, M.,
Melquiond, G.,
Weis, P.:
Trusting computations: a mechanized proof from partial differential
equations to actual program.
Computers & Mathematics with Applications

[18]

Boldo, S.,
Clément, F.,
Filliâtre, J.-C.,
Mayero, M.,
Melquiond, G.,
Weis, P.:
Wave equation numerical resolution: a comprehensive mechanized proof
of a c program.
Journal of Automated Reasoning

[19]

Tekriwal, M.,
Duraisamy, K.,
Jeannin, J.-B.:
A formal proof of the lax equivalence theorem for finite difference
In: Dutle, A.,
Moscato, M.M.,
Titolo, L.,
Muñoz, C.A.,
Perez, I. (eds.)
NASA Formal Methods,
pp. 322–339.

[20]

Immler, F.,
Hölzl, J.:
Numerical analysis of ordinary differential equations in
In: International Conference on Interactive Theorem Proving,
pp. 377–392

[21]

Immler, F.:
A verified ode solver and smale's 14th problem.
Technische Universität München,

[22]

Immler, F.,
Traut, C.:
The flow of odes.
In: International Conference on Interactive Theorem Proving,
pp. 184–199

[23]

Immler, F.,
Traut, C.:
The flow of odes: Formalization of variational equation and
poincaré map.
Journal of Automated Reasoning

[24]

Immler, F.:
Formally verified computation of enclosures of solutions of ordinary
differential equations.
In: Badger, J.M.,
Rozier, K.Y. (eds.)
NASA Formal Methods,
pp. 113–127.

[25]

Mahboubi, A.,
Tassi, E.:
Mathematical components

[26]

Sohrab, H.H.:
Basic Real Analysis.

[27]

The Ratio Test for Sequence Convergence - Mathonline.
(Accessed on 09/21/2021)

[28]

Cohen, C.:
Construction of real algebraic numbers in Coq.
In: Beringer, L.,
Felty, A. (eds.)
ITP - 3rd International Conference on Interactive Theorem Proving -
Princeton, United States

[29]

Cohen, C.,
Mahboubi, A.:
Formal proofs in real algebraic geometry: from ordered fields to
quantifier elimination.
Logical Methods in Computer Science

[30]

Reich, E.:
On the convergence of the classical iterative method of solving linear
simultaneous equations.
The Annals of Mathematical Statistics

[31]

Roux, P.:
Formal proofs of rounding error bounds.
Journal of Automated Reasoning

[32]

Ramananandro, T.,
Mountcastle, P.,
Meister, B.,
Lethin, R.:
A unified coq framework for verifying c programs with floating-point
In: Proceedings of the 5th ACM SIGPLAN Conference on Certified
Programs and Proofs.
CPP 2016,
pp. 15–26.
Association for Computing Machinery,
New York, NY, USA

[33]

Appel, A.W.:
Verified software toolchain.
In: Barthe, G. (ed.)
Programming Languages and Systems,
pp. 1–17.
Berlin, Heidelberg