Since the numerical verification method for solving boundary value problems for elliptic partial differential equations (PDEs) was first developed in 1988, many methods have been devised. In this paper, existing verification methods are reformulated using a convergence theorem for simplified Newton-like methods in the direct product space V-h x V-perpendicular to of a computable finite-dimensional space V-h and its orthogonal complement space V-perpendicular to. Additionally, the Verified Computation for PDEs (VCP) library is provided, which is a software library written in the C++ programming language. The VCP library is introduced as a software library for numerical verification methods of solutions to PDEs. Finally, numerical examples are presented using the reformulated verification methods and VCP library.

<title>Abstract</title>Infinite-dimensional Newton methods can be effectively used to derive numerical proofs of the existence of solutions to partial differential equations (PDEs). In computer-assisted proofs of PDEs, the original problem is transformed into the infinite-dimensional Newton-type fixed point equation <inline-formula><alternatives><tex-math>$$w = - {\mathcal {L } }^{-1} {\mathcal {F } }(\hat{u}) + {\mathcal {L } }^{-1} {\mathcal {G } }(w)$$</tex-math><mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML">
<mml:mrow>
<mml:mi>w</mml:mi>
<mml:mo>=</mml:mo>
<mml:mo>-</mml:mo>
<mml:msup>
<mml:mrow>
<mml:mi>L</mml:mi>
</mml:mrow>
<mml:mrow>
<mml:mo>-</mml:mo>
<mml:mn>1</mml:mn>
</mml:mrow>
</mml:msup>
<mml:mi>F</mml:mi>
<mml:mrow>
<mml:mo>(</mml:mo>
<mml:mover>
<mml:mi>u</mml:mi>
<mml:mo>^</mml:mo>
</mml:mover>
<mml:mo>)</mml:mo>
</mml:mrow>
<mml:mo>+</mml:mo>
<mml:msup>
<mml:mrow>
<mml:mi>L</mml:mi>
</mml:mrow>
<mml:mrow>
<mml:mo>-</mml:mo>
<mml:mn>1</mml:mn>
</mml:mrow>
</mml:msup>
<mml:mi>G</mml:mi>
<mml:mrow>
<mml:mo>(</mml:mo>
<mml:mi>w</mml:mi>
<mml:mo>)</mml:mo>
</mml:mrow>
</mml:mrow>
</mml:math></alternatives></inline-formula>, where <inline-formula><alternatives><tex-math>$${\mathcal {L } }$$</tex-math><mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML">
<mml:mi>L</mml:mi>
</mml:math></alternatives></inline-formula> is a linearized operator, <inline-formula><alternatives><tex-math>$${\mathcal {F } }(\hat{u})$$</tex-math><mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML">
<mml:mrow>
<mml:mi>F</mml:mi>
<mml:mo>(</mml:mo>
<mml:mover>
<mml:mi>u</mml:mi>
<mml:mo>^</mml:mo>
</mml:mover>
<mml:mo>)</mml:mo>
</mml:mrow>
</mml:math></alternatives></inline-formula> is a residual, and <inline-formula><alternatives><tex-math>$${\mathcal {G } }(w)$$</tex-math><mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML">
<mml:mrow>
<mml:mi>G</mml:mi>
<mml:mo>(</mml:mo>
<mml:mi>w</mml:mi>
<mml:mo>)</mml:mo>
</mml:mrow>
</mml:math></alternatives></inline-formula> is a nonlinear term. Therefore, the estimations of <inline-formula><alternatives><tex-math>$$\Vert {\mathcal {L } }^{-1} {\mathcal {F } }(\hat{u}) \Vert $$</tex-math><mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML">
<mml:mrow>
<mml:mrow>
<mml:mo>‖</mml:mo>
</mml:mrow>
<mml:msup>
<mml:mrow>
<mml:mi>L</mml:mi>
</mml:mrow>
<mml:mrow>
<mml:mo>-</mml:mo>
<mml:mn>1</mml:mn>
</mml:mrow>
</mml:msup>
<mml:mi>F</mml:mi>
<mml:mrow>
<mml:mo>(</mml:mo>
<mml:mover>
<mml:mi>u</mml:mi>
<mml:mo>^</mml:mo>
</mml:mover>
<mml:mo>)</mml:mo>
</mml:mrow>
<mml:mrow>
<mml:mo>‖</mml:mo>
</mml:mrow>
</mml:mrow>
</mml:math></alternatives></inline-formula> and <inline-formula><alternatives><tex-math>$$\Vert {\mathcal {L } }^{-1}{\mathcal {G } }(w) \Vert $$</tex-math><mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML">
<mml:mrow>
<mml:mrow>
<mml:mo>‖</mml:mo>
</mml:mrow>
<mml:msup>
<mml:mrow>
<mml:mi>L</mml:mi>
</mml:mrow>
<mml:mrow>
<mml:mo>-</mml:mo>
<mml:mn>1</mml:mn>
</mml:mrow>
</mml:msup>
<mml:mrow>
<mml:mi>G</mml:mi>
<mml:mrow>
<mml:mo>(</mml:mo>
<mml:mi>w</mml:mi>
<mml:mo>)</mml:mo>
</mml:mrow>
<mml:mo>‖</mml:mo>
</mml:mrow>
</mml:mrow>
</mml:math></alternatives></inline-formula> play major roles in the verification procedures . In this paper, using a similar concept to block Gaussian elimination and its corresponding ‘Schur complement’ for matrix problems, we represent the inverse operator <inline-formula><alternatives><tex-math>$${\mathcal {L } }^{-1}$$</tex-math><mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML">
<mml:msup>
<mml:mrow>
<mml:mi>L</mml:mi>
</mml:mrow>
<mml:mrow>
<mml:mo>-</mml:mo>
<mml:mn>1</mml:mn>
</mml:mrow>
</mml:msup>
</mml:math></alternatives></inline-formula> as an infinite-dimensional operator matrix that can be decomposed into two parts: finite-dimensional and infinite-dimensional. This operator matrix yields a new effective realization of the infinite-dimensional Newton method, which enables a more efficient verification procedure compared with existing Nakao’s methods for the solution of elliptic PDEs. We present some numerical examples that confirm the usefulness of the proposed method. Related results obtained from the representation of the operator matrix as <inline-formula><alternatives><tex-math>$${\mathcal {L } }^{-1}$$</tex-math><mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML">
<mml:msup>
<mml:mrow>
<mml:mi>L</mml:mi>
</mml:mrow>
<mml:mrow>
<mml:mo>-</mml:mo>
<mml:mn>1</mml:mn>
</mml:mrow>
</mml:msup>
</mml:math></alternatives></inline-formula> are presented in the “Appendix”.

This paper is concerned with an explicit value of the embedding constant from W-1,W- q(Omega) to L-p(Omega) for a domain Omega subset of R-N (N is an element of N), where 1 &lt;= q &lt;= p &lt;=infinity. We previously proposed a formula for estimating the embedding constant on bounded and unbounded Lipschitz domains by estimating the norm of Stein's extension operator. Although this formula can be applied to a domain Omega that can be divided into a finite number of Lipschitz domains, there was room for improvement in terms of accuracy. In this paper, we report that the accuracy of the embedding constant is significantly improved by restricting Omega to a domain dividable into bounded convex domains.

A fast numerical verification method is proposed for evaluating the accuracy of numerical solutions for symmetric saddle point linear systems whose diagonal blocks of the coefficient matrix are semidefinite matrices. The method is based on results of an algebraic analysis of a block diagonal preconditioning. Some numerical experiments are present to illustrate the usefulness of the method.

This paper focuses on blow-up solutions of ordinary differential equations (ODEs). We present a method for validating blow-up solutions and their blow-up times, which is based on compactifications and the Lyapunov function validation method. The necessary criteria for this construction can be verified using interval arithmetic techniques. Some numerical examples are presented to demonstrate the applicability of our method. (C) 2016 Elsevier B.V. All rights reserved.

In this paper, we propose a verified numerical method for obtaining a sharp inclusion of the best constant for the embedding H-0(1)(Omega) hooked right arrow L-p (Omega) on a bounded convex domain in R-2. We estimate the best constant by computing the corresponding extremal function using a verified numerical computation. Verified numerical inclusions of the best constant on a square domain are presented. (C) 2016 Elsevier B.V. All rights reserved.

In this study, we examine the accurate matrix multiplication in floating-point arithmetic. We demonstrate the error-free transformations of matrix multiplication using high performance basic linear algebra subprograms. These transformations can be applied to accurately compute the product of two matrices using floating-point entries. A key technique for this calculation is error-free splitting in floating-point matrices. In this study, we improve upon our previous method by a posteriori validation using floating-point exception. In the method, we utilize the presence of overflow in a positive manner for detecting whether rounding error occurs. If overflow occurs, the result contains some exceptional values such as +/- and NaN, that is, the method fails by necessity. Otherwise, we can confirm that no rounding error occurs in the process. Therefore, reducing the possibility of overflow is important. The numerical results suggest that the proposed algorithm provides more accurate results compared with the original algorithm. Moreover, for the product of n x n matrices, when n5000, the new algorithm reduces the computing time for error-free transformation by an average of 20 % and up to 30 % compared with the original algorithm. Furthermore, the new algorithm can be used when matrix multiplication is performed using divide-and-conquer methods. Copyright (c) 2016 John Wiley & Sons, Ltd.

This paper is concerned with floating-point filters for a two dimensional orientation problem which is a basic problem in the field of computational geometry. If this problem is only approximately solved by floating-point arithmetic, then an incorrect result may be obtained due to accumulation of rounding errors. A floating-point filter can quickly guarantee the correctness of the computed result if the problem is well-conditioned. In this paper, a simple semi-static floating-point filter which handles floating-point exceptions such as overflow and underflow by only one branch is developed. In addition, an improved fully-static filter is developed.

In this paper, we propose a method for estimating the Sobolev-type embedding constant from W-1,W-q(Omega) to L-p(Omega) on a domain Omega subset of R-n (n = 2,3, ... ) with minimally smooth boundary (also known as a Lipschitz domain), where p is an element of(n/(n - 1), infinity) and q = np/(n + p). We estimate the embedding constant by constructing an extension operator from W-1,W-q(Omega) to W-1,W-q(R-n) and computing its operator norm. We also present some examples of estimating the embedding constant for certain domains.

Recently, new algorithms for accurate matrix multiplication have been developed by the authors. A characteristic of the algorithms is a high dependency on level-3 BIAS routines, which are highly optimized for several architectures. An error-free splitting for floating-point matrices is a key technique in the algorithms. In this paper, an improvement of the error-free splitting is focused on. It is shown by numerical examples that the accuracy of computed results of matrix products can be improved by the modified error-free splitting, compared to that by the previous algorithms. (C) 2015 Elsevier B.V. All rights reserved.

This paper is concerned with real interval arithmetic. We focus on interval matrix multiplication. Well-known algorithms for this purpose require the evaluation of several point matrix products to compute one interval matrix product. In order to save computing time we propose a method that modifies such known algorithm by partially using low-precision floating-point arithmetic. The modified algorithms work without significant loss of tightness of the computed interval matrix product but are about 30% faster than their corresponding original versions. The negligible loss of accuracy is rigorously estimated.

Improved componentwise error bounds for approximate solutions of linear systems are derived in the case where the coefficient of a given linear system is an H-matrix. One of the error bounds presented in this paper proves to be tighter than the existing error bound, which is effective especially for ill-conditioned cases. Numerical experiments are performed to illustrate the effect of the improvements.

This paper is concerned with factorization of symmetric and positive definite matrices which are extremely ill-conditioned. Following the results by Rump (1990), Oishi et al. (2007, 2009) and Ogita (2010), Ogita and Oishi (2012) derived an iterative algorithm for an accurate inverse matrix factorization based on Cholesky factorization for such ill-conditioned matrices. We analyze the behavior of the algorithm in detail and give reasons for convergency by the use of numerical error analysis. Main analysis is that each iteration reduces the condition number of a preconditioned matrix by a factor around the relative rounding error unit until convergence. This behavior is consistent with the numerical results.

This paper proposes a verified numerical method of proving the invertibility of linear elliptic operators. This method also provides a verified norm estimation for the inverse operators. This type of estimation is important for verified computations of solutions to elliptic boundary value problems. The proposed method uses a generalized eigenvalue problem to derive the norm estimation. This method has several advantages. Namely, it can be applied to two types of boundary conditions: the Dirichlet type and the Neumann type. It also provides a way of numerically evaluating lower and upper bounds of target eigenvalues. Numerical examples are presented to show that the proposed method provides effective estimations in most cases.

早稲田大学基幹理工学部に応用数理学科を設立して7年経過した．これは工学と数学を半分づつ学ぶ学科である．3年生に回路理論を工学系の中心的な必修科目として設置した．電気回路基礎から電子回路（アナログ電子回路もディジタル電子回路も含む）までを30回の90分講義で教える科目である．この科目を2012年度と2013年度の2年間担当した．講義の準備に当たって，回路理論を数学的な論理性を保って講義することができるかを考えた．講義開始前の数か月と講義開始後の8か月ぐらいの約1年間でこのことに対する思索（とそれに必要な歴史的な文献調査，外国の教科書の調査，現在の技術動向調査が含まれる）を巡らし，その結果をコロナ社から回路理論として出版した．結果的にこの本はマクスウェルの方程式を公理として仮定し，素子特性は数理モデルとして与えられていると考えて回路理論を数学的な論理性を保つように展開することを志した．教科書とするために，数学的道具立ては制限した部分が多く，また，原稿も半分程度に圧縮したが，回路理論の論理的展開のためにいろいろな講義展開法についての試行を行い，我が国の定石の講義法とかなり異なっている部分も多い．2年目には講義の前半は理論，後半は実験という形で講義を展開した．2年目は回路理論を論理的に捉えるだけでなく，実験により回路の実在をどう捉えるか考えさせた．これらの思索と実践について報告する．

This paper describes the emergence of two characteristic notions, nondegenerate Peano continuum and Cantor set, by the exploration of the essence of the existence of primitive chaos from a topological viewpoint. The primitive chaos is closely related to vital problems in physics itself and leads to chaotic features under natural conditions. The nondegenerate Peano continuum represents an ordinarily observed space, and the existence of a single nondegenerate Peano continuum guarantees the existence of infinite varieties of the primitive chaos leading to the chaos. This result provides an explanation of the reason why we are surrounded by diverse chaotic behaviors. Also, the Cantor set is a general or universal notion different from the special set, the Cantor middle-third set, and the existence of a single Cantor set guarantees infinite varieties of the primitive chaos leading to the chaos. This analogy implies the potential of the Cantor set for the method of new recognizing physical phenomena.

An efficient format and fast algorithms of basic operations for 4-fold working precision are proposed. The proposed format is an unevaluated sum of four double precision numbers, capable of representing at least 203 bits of mantissa. Hence, it is slightly less accurate than quad-double format proposed by Hida et. al. [1], however presented algorithms based on the format are faster than those algorithms. By numerical experiments it is shown that the proposed algorithms are efficient.

This paper presents an algorithm of identifying parameters satisfying a sufficient condition of Plum's Newton-Kantorovich like theorem. Plum's theorem yields a numerical existence test of solutions for nonlinear partial differential equations. The sufficient condition of Plum's theorem is given by the nonemptiness of a region defined by one dimensional nonlinear inequalities. The aim of this paper is to develop a systematic method of constructing an inner inclusion of this region. If $\underline\rho \in {\mathbb R}^+$ is the minimum included in this region, $\underline\rho$ gives the minimum of the error bounds. Moreover, if $\overline\rho \in {\mathbb R}^+$ is the maximum included in this region, then $\overline\rho$ gives the maximum radius of a ball in which the exact solution is unique. In this paper, an algorithm is developed for finding $\rho_e$ and $\rho_u$ such that they belong to this region and become close approximations of $\underline\rho$ and $\overline\rho$, respectively. Finally, to illustrate features of Plum's theorem with our proposed algorithm, some numerical results compared with results by Plum's theorem with the Newton method are presented. In addition to this, Plum's theorem with our algorithm is also compared with Newton-Kantorovich's theorem. One of the most important facts found in this paper is that, for some examples, $\rho_e$ become smaller than error bounds obtained by Newton-Kantorovich's theorem. Moreover, also for these examples, $\rho_u$ become greater than regions indicating uniqueness of the exact solution derived by Newton-Kantorovich's theorem. This implies that Plum's theorem can be seen as a modification of Newton-Kantorovich's theorem.

For Poisson's equation over a polygonal domain of general shape, the solution of which may have a singularity around re-entrant corners, we provide an explicit a priori error estimate for the approximate solution obtained by finite element methods of high degree. The method used herein is a direct extension of the one developed in preceding paper of the second and third listed authors, which provided a new approach to deal with the singularity by using linear finite elements. In the present paper, we also give a detailed discussion of the dependency of the convergence order on solution singularities, mesh sizes and degrees of the finite element method used.

In this paper, a numerical verification method is presented for second-order semilinear elliptic boundary value problems on arbitrary polygonal domains. Based on the Newton-Kantorovich theorem, our method can prove the existence and local uniqueness of the solution in the neighborhood of its approximation. In the treatment of polygonal domains with an arbitrary shape, which gives a singularity of the solution around the re-entrant corner, the computable error estimate of a projection into the finite-dimensional function space plays an essential role. In particular, the lack of smoothness of the solution makes classical error estimates fail on nonconvex domains. By using the Hyper-circle equation, an alternative error estimate of the projection has been proposed. Additionally, a new residual evaluation method based on the mixed finite element method works well. It yields more accurate evaluation than the existing method. The efficiency of our method is shown through illustrative numerical results on several polygonal domains.

This paper is concerned with verification methods for numerical solutions of linear systems. Many methods for the verification require switches of rounding modes defined by the IEEE 754 standard. However, the switches cannot be supported in several computational environments. In such cases, Ogita-Rump-Oishi's method can work on such environments. Recently, Rump developed new error estimates of floating-point summation and dot product. The aim of this paper is to improve Ogita-Rump-Oishi's error estimates by using the error estimates by Rump. In addition, the computational cost of our method is comparable to that of Ogita-Rump-Oishi's method.

This paper is concerned with accurate numerical algorithms for matrix multiplication. Recently, an error-free transformation from a product of two floating-point matrices into an unevaluated sum of floating-point matrices has been developed by the authors. Combining this technique and accurate summation algorithms, new algorithms for accurate matrix multiplication could be investigated. In this paper, it is mentioned that the previous work is not the unique way to achieve an error-free transformation and the constraint of the error-free transformation is clarified. For the application, a new algorithm is developed reducing the number of matrix products compared to the previous algorithm.

The authors recently published a paper on some properties of the solution curves for the last n-1 equations of F(x)+Ax=b where x=[x1,x2,...,xn]T is a variable, F(x):Rn → Rn is a nonlinear function of which the first and second derivatives are strictly positive, A ∈ Rn × n is an Ω-matrix, and b ∈ Rn is a constant vector. In that paper, the authors showed that any solution curve possesses neither maximal points nor inflection points with respect to x1, by making use of a fundamental property of Ω-matrices, which is expressed in the form of inequality. However, the proof was a little unclear as shown in Introduction. The objective of this paper is to give an explicit proof for the property of Ω-matrices, which makes the author's previous result more rigorous.

The finite element method (FEM) is applied to bound leading eigenvalues of the Laplace operator over polygonal domains. Compared with classical numerical methods, most of which can only give concrete eigenvalue bounds over special domains of symmetry, our proposed algorithm can provide concrete eigenvalue bounds for domains of arbitrary shape, even when the eigenfunction has a singularity. The problem of eigenvalue estimation is solved in two steps. First, we construct a computable a priori error estimation for the FEM solution of Poisson's problem, which holds even for nonconvex domains with reentrant corners. Second, new computable lower bounds are developed for the eigenvalues. Because the interval arithmetic is implemented throughout the computation, the desired eigenvalue bounds are expected to be mathematically correct. We illustrate several computation examples, such as the cases of an L-shaped domain and a crack domain, to demonstrate the efficiency and flexibility of the proposed method.

We consider an explicit estimation for error constants from two basic constant interpolations on triangular finite elements. The problem of estimating the interpolation constants is related to the eigenvalue problems of the Laplacian with certain boundary conditions. By adopting the Lehmann-Goerisch theorem and finite element spaces with a variable mesh size and polynomial degree, we succeed in bounding the leading eigenvalues of the Laplacian and the error constants with high precision. An online demo for the constant estimation is also available at http://www.xfliu.org/onlinelab/.

Since the chaos was discovered, it has been recognized that we are surrounded by diverse chaotic behaviors. The purpose of this study is to reconsider the implication of this fact through the notion of a primitive chaos. Under natural conditions, each primitive chaos leads to apparent chaotic features, such as the existence of a nonperiodic orbit, the existence of the periodic point whose prime period is n for any n is an element of N, the existence of a dense orbit, the density of periodic points, sensitive dependence on initial conditions, and topological transitivity, while infinite varieties of the primitive chaos are guaranteed by a nondegenerate Peano continuum.

This paper concerns a robust algorithm for the 2D orientation problem which is one of the basic tasks in computational geometry. Recently, a fast and accurate floating-point summation algorithm is investigated by Rump, Ogita and Oishi in [S.M. Rump, T. Ogita, S. Oishi, Accurate floating-point summation. Part I: Faithful rounding, SIAM J. Sci. Comput. 31 (1) (2008) 189-2241, in which a new kind of an error-free transformation of floating-point numbers is used. Based on it, a new algorithm of error-free determinant transformation for the 2D orientation problem is proposed, which gives a correct result. Numerical results are presented for illustrating that the proposed algorithm has some advantage over preceding algorithms in terms of measured computing time. (C) 2012 Elsevier Inc. All rights reserved.

This paper is concerned with accurate matrix multiplication in floating-point arithmetic. Recently, an accurate summation algorithm was developed by Rump et al. (SIAM J Sci Comput 31(1):189-224, 2008). The key technique of their method is a fast error-free splitting of floating-point numbers. Using this technique, we first develop an error-free transformation of a product of two floating-point matrices into a sum of floating-point matrices. Next, we partially apply this error-free transformation and develop an algorithm which aims to output an accurate approximation of the matrix product. In addition, an a priori error estimate is given. It is a characteristic of the proposed method that in terms of computation as well as in terms of memory consumption, the dominant part of our algorithm is constituted by ordinary floating-point matrix multiplications. The routine for matrix multiplication is highly optimized using BLAS, so that our algorithms show a good computational performance. Although our algorithms require a significant amount of working memory, they are significantly faster than 'gemmx' in XBLAS when all sizes of matrices are large enough to realize nearly peak performance of 'gemm'. Numerical examples illustrate the efficiency of the proposed method.

We discuss several methods for real interval matrix multiplication. First, earlier studies of fast algorithms for interval matrix multiplication are introduced: naive interval arithmetic, interval arithmetic by midpoint-radius form by Oishi-Rump and its fast variant by Ogita-Oishi. Next, three new and fast algorithms are developed. The proposed algorithms require one, two or three matrix products, respectively. The point is that our algorithms quickly predict which terms become dominant radii in interval computations. We propose a hybrid method to predict which algorithm is suitable for optimizing performance and width of the result. Numerical examples are presented to show the efficiency of the proposed algorithms. (C) 2011 Elsevier B.V. All rights reserved.

This paper is concerned with an accurate computing for matrix multiplication. We show that matrix multiplication can be transformed into a summation of floating-point matrices by mainly using level 3 operations in BLAS. We call it 'error-free transformation of matrix multiplication'. By combining this error-free transformation and accurate summation algorithms for floating-point numbers, we can obtain an accurate result of matrix multiplication. Numerical examples are presented to illustrate the efficiency of the proposed algorithm.

In this paper, an algorithm for an accurate matrix factorization based on Cholesky factorization for extremely ill-conditioned matrices is proposed. The Cholesky factorization is widely used for solving a system of linear equations whose coefficient matrix is symmetric and positive definite. However, it sometimes breaks down by the presence of an imaginary root due to the accumulation of rounding errors, even if the matrix is actually positive definite. To overcome this, a completely stable algorithm named inverse Cholesky factorization is investigated, which never breaks down as long as the matrix is symmetric and positive definite. The proposed algorithm consists of standard numerical algorithms and an accurate algorithm for dot products. Moreover, it is shown that the algorithm can also verify the positive definiteness of a given real symmetric matrix. Numerical results are presented for illustrating the performance of the proposed algorithms.

In this paper we study the generation of an ill-conditioned integer matrix A=[aij] with |aij|≤µ for some given constant µ. Let n be the order of A. We first give some upper bounds of the condition number of A in terms of n and µ. We next propose new methods to generate extremely ill-conditioned integer matrices. These methods are superior to the well-known method by Rump in some respects, namely, the former has a simple algorithm to generate a larger variety of ill-conditioned matrices. In particular we propose a method to generate ill-conditioned matrices with a choice of desirable singular value distributions as benchmark matrices.

This paper is concerned with the tight enclosure of matrix multiplication AB for two floating-point matrices A and B. The aim of this paper is to compute component-wise upper and lower bounds of the exact result C of the matrix multiplication AB by floating-point arithmetic. Namely, an interval matrix enclosing C is obtained. In this paper, new algorithms for enclosing C are proposed. The proposed algorithms are designed to mainly exploit the level 3 operations in BLAS. Although the proposed algorithms take around twice as much costs as a standard algorithm promoted by Oishi and Rump, the accuracy of the result by the proposed algorithms is better than that of the standard algorithm. At the end of this paper, we present numerical examples showing the efficiency of the proposed algorithms. Copyright (C) 2010 John Wiley & Sons, Ltd.

Several methods have been proposed to calculate a rigorous error bound of an approximate solution of a linear system by floating-point arithmetic. These methods are called &apos;verification methods&apos;. Applicable range of these methods are different. It depends mainly on the condition number and the dimension of the coefficient matrix whether such methods succeed to work or not. In general, however, the condition number is not known in advance. If the dimension or the condition number is large to some extent, then Oishi-Rump&apos;s method, which is known as the fastest verification method for this purpose, may fail. There are more robust verification methods whose computational cost is larger than the Oishi-Rump&apos;s one. It is not so efficient to apply such robust methods to well-conditioned problems. The aim of this paper is to choose a suitable verification method whose computational cost is minimum to succeed. First in this paper, four fast verification methods for linear systems are briefly reviewed. Next, a compromise method between Oishi-Rump&apos;s and Ogita-Oishi&apos;s one is developed. Then, an algorithm which automatically and efficiently chooses an appropriate verification method from five verification methods is proposed. The proposed algorithm does as much work as necessary to calculate error bounds of approximate solutions of linear systems. Finally, numerical results are presented.

This paper is concerned with interval arithmetic, especially, an enclosure of a matrix product is focused on. Using level 3 operations of matrix computations, an algorithm outputting a tight enclosure for matrix multiplication is proposed. Most of the algorithms for this purpose require switches of rounding modes defined in the IEEE standard 754. However some programing enviroments have not supported them. Our proposed method demands only rounding-to-nearest mode, so that it is very portable.

Present authors have presented with Takayuki Kubo at University of Tsukuba a method of a computer assisted proof for the existence and uniqueness of solutions to two-point boundary value problems of nonlinear ordinary differential equations in the paper submitted for NOLTA, IEICE. This method uses piecewise linear finite element base functions and sometimes requires fine mesh. To overcome this difficulty, in this paper, an improved method is presented for the norm estimation of the residual to the operator equation. In this refined formulation, piecewise quadratic finite element base functions are used. A kind of the residual technique works sophisticatedly well. It is stated that the estimation of the residual can be expected smaller than that of the previous method. Finally, four examples are presented. Each result demonstrates that a remarkable improvement is achieved in accuracy of the guaranteed error estimation.

It is well known that it is an ill-posed problem to decide whether a function has a multiple root. For example, an arbitrarily small perturbation of a real polynomial may change a double real root into two distinct real or complex roots. In this paper we describe a computational method for the verified computation of a complex disc to contain exactly k roots of a univariate nonlinear function. The function may be given by some program. Computational results using INTLAB, the Matlab toolbox for reliable computing, demonstrate properties and limits of the method.

Given a vector pi of floating-point numbers with exact sum s, we present a new algorithm with the following property: Either the result is a faithful rounding of s, or otherwise the result has a relative error not larger than epsKcond(∑pi) for K to be specified. The statements are also true in the presence of underflow, the computing time does not depend on the exponent range, and no extra memory is required. Our algorithm is fast in terms of measured computing time because it allows good instruction-level parallelism. A special version for K=2, i.e., quadruple precision is also presented. Computational results show that this algorithm is more accurate and faster than competitors such as XBLAS.

A fast verified automatic integration algorithm is proposed for calculating univariate integrals over finite intervals. This algorithm is based on the double exponential formula proposed by Takahasi and Mori. The double exponential formula uses a certain trapezoidal rule. This trapezoidal rule is determined by fixing two parameters, the width h of a subdivision of a finite interval and the number n of subdivision points of this subdivision. A theorem is presented for calculating h and n as a function of a given tolerance of the verified numerical integration of a definite integral. An efficient a priori method is also proposed for evaluating function calculation errors including rounding errors of floating point calculations. Combining these, a fast algorithm is proposed for verified automatic integration. Numerical examples are presented for illustrating effectiveness of the proposed algorithm.

In this paper, a numerical method is presented for verifying the existence and uniqueness of solutions to two-point boundary value problems of nonlinear ordinary differential equations. Taking into account every error of numerical computations such as the discretization error and the rounding error, this method also provides mathematically guaranteed error bounds between approximations obtained by numerical computations and the exact solution whose existence is proven by the numerical existence theorem, which is based on the Newton-Kantorovich theorem. Finally, illustrative numerical results are presented for showing the usefulness of the method.

A verified automatic integration algorithm is proposed for calculating contour integration over complex field using numerical computations. The proposed algorithm is based on trapezoidal rule for angle. The error analysis of the method have been presented by several authors, however, these investigations are done basically for examining the rates of convergence, and several constants in these error formula were left unevaluated. In order to construct verified numerical integrator using the algorithm, the error formula is presented. To construct efficient verified numerical integrator, an efficient a priori method of evaluating function calculation errors is adopted. Combining these, a verified automatic integration algorithm is proposed. Numerical examples are presented for illustrating effectiveness of the proposed algorithm.

This paper treats a linear equation
Av = b,
where A is an element of F(nxn) and b is an element of F(n). Here, F is a set of floating point numbers. Let u be the unit round-off of the working precision and kappa(A) = parallel to A parallel to(infinity)parallel to A(-1)parallel to(infinity) be the condition number of the problem. In this paper, ill-conditioned problems with
1 &lt; u kappa(A) &lt; infinity
are considered and an iterative refinement algorithm for the problems is proposed. In this paper, the forward and backward stability will be shown for this iterative refinement algorithm.

This paper is concerned with a robust geometric predicate for the 2D orientation problem. Recently, a fast and accurate floating-point summation algorithm is investigated by Rump, Ogita and Oishi, which provably outputs a result faithfully rounded from the exact value of the summation of floating-point numbers. We optimize their algorithm for applying it to the 2D orientation problem which requires only a correct sign of a determinant of a 3 x 3 matrix. Numerical results illustrate that our algorithm works fairly faster than the state-of-the-art algorithm in various cases.

This paper aims to survey fast methods of verifying the accuracy of a numerical solution of a linear system. For the last decade, a number of fast verification algorithms have been proposed to obtain an error bound of a numerical solution of a dense or sparse linear system. Such fast algorithms rely on the verified numerical computation using floating-point arithmetic defined by IEEE standard 754. Some fast verification methods for dense and sparse linear systems are reviewed together with corresponding numerical results to show the practical use and efficiency of the verified numerical computation as much as possible.

This paper concerns with the following linear programming problem: \[ \mbox{Maximize } c^tx, \mbox{ subject to } Ax \leqq b \mbox{ and } x\geqq 0, \] where $A \in \F^{m\times n}$, $b \in \F^m$ and $c, x \in \F^n$. Here, $\F$ is a set of floating point numbers. The aim of this paper is to propose a numerical method of including an optimum point of this linear programming problem provided that a good approximation of an optimum point is given. The proposed method is base on Kantorovich's theorem and the continuous Newton method. Kantorovich's theorem is used for proving the existence of a solution for complimentarity equation and the continuous Newton method is used to prove feasibility of that solution. Numerical examples show that a computational cost to include optimum point is about 4 times than that for getting an approximate optimum solution.

Parallel algorithms for accurate summation and dot product are proposed, They are parallelized versions of fast and accurate algorithms of calculating sum and dot product using error-free transformations which are recently proposed by Ogita et al. [T. Ogita, S.M. Rump, S. Oishi, Accurate sum and dot product, SIAM J. Sci. Comput. 26 (6) (2005) 1955-1988]. They have shown their algorithms are fast in terms of measured computing time. However, due to the strong data dependence in the process of their algorithms, it is difficult to parallelize them. Similarly to their algorithms, the proposed parallel algorithms in this paper are designed to achieve the results as if computed in K-fold working precision with keeping the fastness of their algorithms. Numerical results are presented showing the performance of the proposed parallel algorithm of calculating dot product. (C) 2008 Elsevier B.V. All rights reserved.

筆者は1976 年の卒論より研究をスタートしました。すでに32 年間研究に携わってきたことになります。筆者が精度保証付き数値計算の研究に移ったのは1990年のことです。以来本分野で研究を続けてきました。精度保証付き数値計算の研究の研究に移ったのは筆者なりの必然性があります。1990年当時の精度保証付き数値計算の研究は実用的ではないと考えられていたような気がします。実際、数百次元の連立一次方程式の精度保証が精一杯の感じでありました。現在では特殊な構造を持つ方程式であれば一億次元の連立一次方程式でも精度保証できるようになり、精度保証付き数値計算は実用の段階に至っていると思っています。筆者の研究がこのようなブレークスルーに貢献できたと考えておりますが、本稿ではこのような精度保証付き数値計算の研究の発展と筆者の研究の個人史の交錯を描かせていただきました。

In Part II of this paper we first refine the analysis of error-free vector transformations presented in Part I. Based on that we present an algorithm for calculating the rounded-to-nearest result of s := Sigma pi for a given vector of floating-point numbers pi, as well as algorithms for directed rounding. A special algorithm for computing the sign of s is given, also working for huge dimensions. Assume a floating-point working precision with relative rounding error unit eps. We define and investigate a K-fold faithful rounding of a real number r. Basically the result is stored in a vector Res(nu) of K nonoverlapping floating-point numbers such that Sigma Res(nu) approximates r with relative accuracy eps(K), and replacing Res(K) by its floating-point neighbors in Sigma Res(nu) forms a lower and upper bound for r. For a given vector of floating-point numbers with exact sum s, we present an algorithm for calculating a K-fold faithful rounding of s using solely the working precision. Furthermore, an algorithm for calculating a faithfully rounded result of the sum of a vector of huge dimension is presented. Our algorithms are fast in terms of measured computing time because they allow good instruction-level parallelism, they neither require special operations such as access to mantissa or exponent, they contain no branch in the inner loop, nor do they require some extra precision. The only operations used are standard floating-point addition, subtraction, and multiplication in one working precision, for example, double precision. Certain constants used in the algorithms are proved to be optimal.

Consider the boundary value problem Lu equivalent to -(pu ')' + qu ' + ru = f, a &lt;= x &lt;= b, u(a) = u(b) = 0. Let H(v)A(v)U = f and (A) over cap U-v = (f) over cap be its finite difference equations and piecewise linear finite element equations on partitions Delta(v) : a = x(0)(v) &lt; x(1)(v) &lt; ... &lt; x(nv+1)(v) = b, v = 1, 2.... with h(i)(v) = x(i)(v) - x(i-1)(v), h(v) = max(i)h(i)(v) -&gt; 0 as v -&gt; infinity, where H-v, are n(v) x n(v) diagonal matrices and A(v) as well as (A) over cap (v) are n(v) x n(v) tridiagonal. It is shown that the following three conditions are equivalent: (i) The boundary value problem has a unique solution u is an element of C-2[a, b]. (ii) For sufficiently large v &gt;= v(0), the inverse A(v)(-1) = (g(ij)(v)) exists and vertical bar g(ij)(v)vertical bar &lt;= M, for all(i,j) with a constant M &gt; 0 independent of h(v). (iii) For sufficiently large v &gt;= (v) over cap (0), (A(v)) over cap (-1) = ((g(ij)) over cap)(v) exists and vertical bar(g(ij)) over cap (v)vertical bar &lt;= (M) over cap, for all(i,j) with a constant (M) over cap &gt; 0 independent of h(v). It is also shown by a numerical example that the finite difference method with uniform nodes x(i+1) = x(i) + h, 0 &lt;= i &lt;= n, h = (b - a)/(n + 1) applied to the boundary value problem with no solution gives a ghost solution for every n.

Given a vector of floating-point numbers with exact sum s, we present an algorithm for calculating a faithful rounding of s, i.e., the result is one of the immediate floating-point neighbors of s. If the sum s is a floating-point number, we prove that this is the result of our algorithm. The algorithm adapts to the condition number of the sum, i.e., it is fast for mildly conditioned sums with slowly increasing computing time proportional to the logarithm of the condition number. All statements are also true in the presence of underflow. The algorithm does not depend on the exponent range. Our algorithm is fast in terms of measured computing time because it allows good instruction-level parallelism, it neither requires special operations such as access to mantissa or exponent, it contains no branch in the inner loop, nor does it require some extra precision: The only operations used are standard floating-point addition, subtraction, and multiplication in one working precision, for example, double precision. Certain constants used in the algorithm are proved to be optimal.

Given a vector of floating-point numbers with exact sum s, we present an algorithm for calculating a faithful rounding of s, i.e., the result is one of the immediate floating-point neighbors of s. If the sum s is a floating-point number, we prove that this is the result of our algorithm. The algorithm adapts to the condition number of the sum, i.e., it is fast for mildly conditioned sums with slowly increasing computing time proportional to the logarithm of the condition number. All statements are also true in the presence of underflow. The algorithm does not depend on the exponent range. Our algorithm is fast in terms of measured computing time because it allows good instruction-level parallelism, it neither requires special operations such as access to mantissa or exponent, it contains no branch in the inner loop, nor does it require some extra precision: The only operations used are standard floating-point addition, subtraction, and multiplication in one working precision, for example, double precision. Certain constants used in the algorithm are proved to be optimal.

Consider the boundary value problem Lu equivalent to -(pu ')' + qu ' + ru = f, a &lt;= x &lt;= b, u(a) = u(b) = 0. Let H(v)A(v)U = f and (A) over cap U-v = (f) over cap be its finite difference equations and piecewise linear finite element equations on partitions Delta(v) : a = x(0)(v) &lt; x(1)(v) &lt; ... &lt; x(nv+1)(v) = b, v = 1, 2.... with h(i)(v) = x(i)(v) - x(i-1)(v), h(v) = max(i)h(i)(v) -&gt; 0 as v -&gt; infinity, where H-v, are n(v) x n(v) diagonal matrices and A(v) as well as (A) over cap (v) are n(v) x n(v) tridiagonal. It is shown that the following three conditions are equivalent: (i) The boundary value problem has a unique solution u is an element of C-2[a, b]. (ii) For sufficiently large v &gt;= v(0), the inverse A(v)(-1) = (g(ij)(v)) exists and vertical bar g(ij)(v)vertical bar &lt;= M, for all(i,j) with a constant M &gt; 0 independent of h(v). (iii) For sufficiently large v &gt;= (v) over cap (0), (A(v)) over cap (-1) = ((g(ij)) over cap)(v) exists and vertical bar(g(ij)) over cap (v)vertical bar &lt;= (M) over cap, for all(i,j) with a constant (M) over cap &gt; 0 independent of h(v). It is also shown by a numerical example that the finite difference method with uniform nodes x(i+1) = x(i) + h, 0 &lt;= i &lt;= n, h = (b - a)/(n + 1) applied to the boundary value problem with no solution gives a ghost solution for every n.

In this paper, the problem of inverting regular matrices with arbitrarily large condition number is treated in double precision defined by IEEE 754 floating point standard. In about 1984, Rump derived a method for inverting arbitrarily ill-conditioned matrices. The method requires the possibility to calculate a dot product in higher precision. Rump's method is of theoretical interest. Rump made it clear that inverting an arbitrarily ill-conditioned matrix in single or double precision does not produce meaningless numbers, but contains a lot of information in it. Rump's method uses such inverses as preconditioners. Numerical experiments exhibit that Rump's method converges rapidly for various matrices with large condition numbers. Why Rump's method is so efficient for inverting arbitrarily ill-conditioned matrices is a little mysterious. Thus, to prove its convergence is an interesting problem in numerical error analysis. In this article, a convergence theorem is presented for a variant of Rump's method. (C) 2006 Elsevier B.V. All rights reserved.

Recent development of Java's optimization techniques makes Java one of the most useful programming languages for numerical computations. This paper proposes a numerical method of obtaining verified approximate solutions of linear systems. Usual methods for verified computations use switches of rounding modes defined in IEEE standard 754. However, such switches of rounding modes have not been supported in Java. This method avoids using directed rounding, so that it is implementable on a wide range of programming languages including Java. Numerical experiments using Java illustrate that the method can give a very accurate error bound for an approximate solution of a linear system with almost same computational cost as that for calculating an approximate inverse by the Gaussian elimination. (c) 2005 Elsevier B.V. All rights reserved.

This paper is concerned with semilinear tow-point boundary value problems of the form -(p(x)u')'+ f (x, u) = 0, a &lt;= x &lt;= b, alpha(0)u(a) - alpha(1)u'(a) = alpha, beta(0)u(b) + beta(1)u'(b) = beta, alpha(i) &gt;= 0, beta(i) &gt;= 0, i = 0, 1, alpha(0)+alpha(1) &gt; 0, beta(0)+beta(1) &gt; 0, alpha(0)+beta(0) &gt; 0. Under the assumption inf f(u) &gt; - lambda(1), where lambda(1) is the smallest eigenvalue of Lu = -(pu')' with the boundary conditions, unique existence theorems of solution for the continuous problem and a discretized. system with not necessarily uniform nodes are given as well as error estimates. The results generalize three theorems of Lees for u" = f (x, u), 0 &lt;= x &lt;= 1, u(0) = alpha, u(1) = beta.

This paper is concerned with interval matrix multiplication.New algorithms are proposed to calculate an inclusion of the product of interval matrices using rounding mode controlled computation. Thecomputational cost of the proposed algorithms is almost the same as that for calculating an inclusion of the product of point matrices.Numerical results are presented to illustrate that the new algorithms are much faster than the conventional algorithms and that the guaranteed accuracies obtained by the proposed algorithms are comparable to those of the conventional algorithms. © Springer 2005.

A fast verification algorithm of calculating guaranteed error bounds for all approximate eigenvalues of a real symmetric matrix is proposed. In the proposed algorithm, Rump's and Wilkinson's bounds are combined. By introducing Wilkinson's bound, it is possible to improve the error bound obtained by the verification algorithm based on Rump's bound with a small additional cost. Finally, this paper includes some numerical examples to show the efficiency of the proposed algorithm.

Algorithms for summation and dot product of floating-point numbers are presented which are fast in terms of measured computing time. We show that the computed results are as accurate as if computed in twice or K-fold working precision, K &gt;= 3. For twice the working precision our algorithms for summation and dot product are some 40% faster than the corresponding XBLAS routines while sharing similar error estimates. Our algorithms are widely applicable because they require only addition, subtraction, and multiplication of floating-point numbers in the same working precision as the given data. Higher precision is unnecessary, algorithms are straight loops without branch, and no access to mantissa or exponent is necessary.

This paper proposes the method of numerical verification of solutions of a periodic integral equation with a logarithmic singular kernel, which is typically found in some two-dimensional potential problems. The verification method utilizes a property of the singular integral for trigonometric polynomials, the periodic Sobolev space and Schauder's fixed point theorem.

This paper describes numerical verification of a double turning point of a nonlinear system using an extended system. To verify the existence of a double turning point, we need to prove that one of the solutions of the extended system corresponds to the double turning point. For that, we propose an extended system with an additional condition. As an example, for a finite dimensional problem, we verify the existence and local uniqueness of a double turning point numerically using the extended system and a verification method based on the Banach fixed point theorem.

As interval analysis-based reliable computations find wider application, more software is becoming available. Simultaneously. the applications for which this software is designed are becoming more diverse. Because of this, the software itself takes diverse forms, ranging from libraries for application development to fully interactive systems. The target applications range from fairly general to specialized.
Here, we describe the design of four freely available software systems providing validated computations. Oishi provides Slab, a complete, high-performance system for validated linear algebra whose user interface mimics both Matlab's M-files and a large subset of Matlab's command-line functions. In contrast, CADNA (Fabien Rico) is a C++ library designed to give developers of embedded systems access to validated numeric computations. Addressing global constrained optimization and validated solution of nonlinear algebraic systems, Kearfott's GlobSol focuses on providing the most practical such system possible without specifying non-general problem structure, Kearfott's system has a Fortran-90 interface. Finally, Neher provides a mathematically sound stand-alone package ACETAF with an intuitive graphical user interface for computing complex Taylor coefficients and their bounds, radii of convergence, etc.
Overviews of each package's capabilities, use, and instructions for obtaining and installing appear.

This paper is concerned with the problem of verifying the accuracy of approximate solutions of systems of linear equations. Recently, fast algorithms for calculating guaranteed error bounds of computed solutions of system's of linear equations have been proposed using the rounding mode controlled verification method and the residual iterative verification method. In this paper, a new verification method for systems of linear equations is proposed. Using this verification method, componentwise verified error bounds of approximate solutions of systems of linear equations can be calculated. Numerical results are presented to illustrate that it is possible to get very sharp error bounds of computed solutions of systems of linear equations whose coefficient matrices are symmetric and positive definite.

In this paper. we are concerned with a matrix equation
Ax = b
where A is an a x n real matrix and x and b are n-vectors. Assume that an approximate solution (x) over tilde is given together with an approximate LU decomposition. We will present fast algorithms for proving nonsingularity of A and for calculating rigorous error bounds for parallel toA(-1)b-(x) over tilde parallel to(infinity), The emphasis is on rigour of the bounds. The purpose of this paper is to propose different algorithms, the fastest with 2/3n(3) flops computational cost for the verification step, the same as for the LU decomposition. The presented algorithms exclusively use library routines for LU decomposition and for all other matrix and vector operations.

This paper is concerned with the validation of simple turning points of two-point boundary value problems of nonlinear ordinary differential equations. Usually it is hard to validate approximate solutions of turning points numerically because of it's singularity. In this paper, it is pointed out that applying the infinite dimensional Krawcyzk-based interval validation method to enlarged system, the existence of simple turning points can be verified. Taking an example, the result of validation is also presented.

Synchronous distributed timing clocks are the basic building blocks in digital communication systems. Conventional systems mainly employ a tree-like network of cascaded timing clocks for synchronous clocking. On the other hand, decentralized synchronous networks of timing clocks, which have been proposed from a very early stage of the digital communication, are gaining attention in the consumer communication networks and also recently in large, high-performance digital systems (such as multiprocessors) clocking. In this paper, we present a theoretical study of synchronous networks of timing clocks consisting of locally connected second order phase-locked loops (PLLs). We find a close connection between the stability properties of the first and second order networks. The particular examples of one way and two way nearest neighbor coupling, with a lag-lead filter and a triangular phase detector (PD) are analyzed in detail. Both the synchronized in-phase solution and the wave-like ''mode-lock'' solution are examined. A criterion is found for the stability of the one-way coupled network while the two-way coupled network is found to be always stable.

We analyze the collective behavior of a set of coupled damped driven pendula with finite (large) inertia, and show that the synchronization of the oscillators exhibits a first order phase transition synchronization onset, substantially different from the second order transition obtained in the case of no inertia. There is hysteresis between two macroscopic states, a weakly and a strongly coherent synchronized state, depending on the coupling and the initial state of the oscillators. A self-consistent theory is shown to determine these cooperative phenomena and to predict the observed numerical data in specific examples.

We analyze a large system of nonlinear phase oscillators with sinusoidal nonlinearity, uniformly distributed natural frequen cies and global all-to-all coupling, which is an extension of Kuramoto's model to second-order systems. For small coupling, the system evolves to an incoherent state with the phases of all the oscillators distributed uniformly. As the coupling is increased, the system exhibits a discontinuous transition to the coherently synchronized state at a pinning threshold of the coupling strength, or to a partially synchronized oscillation coherent state at a certain threshold below the pinning threshold. if the coupling is decreased from a strong coupling with all the oscillators synchronized coherently, this coherence can persist until the depinning threshold which is less than the pinning threshold, resulting in hysteretic synchrony depending on the initial configuration of the oscillators. We obtain analytically both the pinning and depinning threshold and also explain the discontinuous transition at the thresholds for the underdamped case in the large system size limit. Numerical exploration shows the oscillatory partially coherent state bifurcates at the depinning threshold and also suggests that this state persists independent of the system size. The system studied here provides a simple model for collective behaviour in damped driven high-dimensional Hamiltonian systems which can explain the synchronous firing of certain fireflies or neural oscillators with frequency adaptation and may also be applicable to interconnected power systems.

Dynamical properties such as lock-in or out-of-lock condition of mutually coupled phase-locked loops (PLL's) are problems of practical interest, The present paper describes a study of such dynamical properties for mutually coupled PLL's incorporating lag filters and triangular phase detectors, The fourth-order ordinary differential equation (ODE) governing the mutually coupled PLL's is reduced to the equivalent third-order ODE due to the symmetry, where the system is analyzed in the context of nonlinear dynamical system theory, An understanding as to how and when lock-in can be obtained or out-of-lock behavior persists, is provided by the geometric structure of the invariant manifolds generated in the vector field from the third-order ODE. In addition, a connection to the recently developed theory on chaos and bifurcations from degenerated homoclinic points is also found to exist. The two-parameter diagrams of the one-homoclinic orbit are obtained by graphical solution of a set of nonlinear (finite dimensional) equations. Their graphical results useful in determining whether the system undergoes lock-in or continues out-of-lock behavior, are verified by numerical simulations.