# From Penrose Equations to Zhang Neural Network, Getz–Marsden Dynamic System, and DDD (Direct Derivative Dynamics) Using Substitution Technique

Authors: Dongqing Wu; Yunong Zhang
Journal: Discrete Dynamics in Nature and Society
 (2021)

---

## Abstract

The time-variant matrix inversion (TVMI) problem solving is the hotspot of current research because of its frequent appearance and application in scientific research and industrial production. The generalized inverse problem of singular square matrix and nonsquare matrix can be related to Penrose equations (PEs). The PEs implicitly define the generalized inverse of a known matrix, which is of fundamental theoretical significance. Therefore, the in-depth study of PEs might enlighten problem solving of TVMI in a foreseeable way. For the first time, we construct three different matrix error-monitoring functions based on PEs and propose the corresponding models for TVMI problem solving by using the substitution technique and ZNN design formula. In order to facilitate computer simulation, the obtained continuous-time models are discretized by using ZTD (Zhang time discretization) formulas. Furthermore, the feasibility and effectiveness of the novel Zhang neural network (ZNN) multiple-multiplication model for matrix inverse (ZMMMI) and the PEs-based Getz–Marsden dynamic system (PGMDS) model in solving the problem of TVMI are investigated and shown via theoretical derivation and computer simulation. Computer experiment results also illustrate that the direct derivative dynamics model for TVMI is less effective and feasible.


---

## Body

## 1. Introduction

During the past decades, scientists and engineers have encountered linear matrix equation problems, i.e., Lyapunov equation, Sylvester equation, and the variational problem, again and again in various scenarios. The matrix inversion problem is one of the most prominent subproblems in the linear matrix equation problems. Among the problems encountered in a variety of optimization problems, the fundamental one is the solution of matrix inversion, such as signal processing [1], biomedical prediction [2], image reconstruction [3], nonlinear optimization [4], and robot inverse kinematics [5–8]. Generally speaking, the matrix inversion problem can be formulated as AX=I, where A∈ℝn×n is a known constant matrix, X∈ℝn×n is the unknown matrix to be computed, and I∈ℝn×n is an identity matrix. Over the years, efforts were directed towards computational issues of time-variant matrix inversion (TVMI) and a wealth of algorithms were proposed and applied to solving this problem [9–11]. For example, Yeung and Kumbi [9] developed an inversion method based on the multidimensional discrete Fourier transform of matrix sequences and applied it to an electron amplifier. In [10], Benner et al. adopted Gauss–Huard algorithm to solve TVMI problem. In [11], Xiao et al. designed a complex-valued nonlinear recurrent neural network to solve time-varying complex matrix inversion. Zhang et al. [12] discussed the solution of TVMI problem based on direct derivative dynamics (DDD) and proposed a zero-stable, consistent, and convergent four-step model.According to the classic solution of matrix inversion, the time complexity of arithmetic operations is proportional toOn3, where n represents the dimension of the square matrix [13]. Evidently, it is an unbearable cost for solving the TVMI problem. Therefore, it is an urgent need for effective methods in TVMI problem solving.In recent years, with the rise of the artificial neural network, the recurrent neural network (RNN) has been thought and manifested to be an effective option for TVMI problem solving. Methods based on RNN have been proved to be promising for RNN’s nature of high-speed parallel-processing and convenience of hardware implementation [14–18]. To solve the TVMI problem, Zhang et al. firstly proposed a new class of RNN, Zhang neural network, abbreviated as ZNN [19]. The essence of ZNN model is to construct a carefully selected error-monitoring function, termed Zhang function (ZF). The ZF can be negative, zero, positive, bounded, or even lower-unbounded. The traditional error-monitoring functions are often some norm-based positive-definite energy functions, which are not so flexible as ZFs. Moreover, since the time-derivative information of the time-variant matrix is fully utilized by ZFs, the resultant ZNN model can decrease the lagging errors which are almost unavoidable in the traditional methods. In theoretical and practical researches, the ZNN models claim their natural advantages in convergence and high accuracy. Inspired by this idea, the ZFs have been found to speed up and consolidate the development of various ZNN models [20–26].With the deepening of the research on solving TVMI, designing a new ZF through specific classic equations, deriving new solution models, and discussing the convergence and accuracy of the models have become the key to the solution of the extended TVMI problem.As a generalization of matrix inverse, the Moore–Penrose generalized inverse has been widely investigated [27]. It is applied to finding the least norm square solution of nonuniform linear equations and makes the form of solution simple. The Moore–Penrose generalized inverse of a matrix is unique in both real and complex fields, which can be obtained theoretically by singular value decomposition (SVD) algorithm [27]. Although the SVD algorithm provides a direction for TVMI problem solving, it is deeply troubled by the high time complexity of the algorithm. For example, according to the actual number of multiplication operations, SVD algorithm needs 2mn2+4n3 operations to find the Moore–Penrose generalized inverse of an m×n matrix [27]. Therefore, based on the mathematical definition of inverse matrix, there have been many attempts to solve the TVMI problem [28–30].In view of the fundamental theoretical value of Penrose equations (PEs) in the definition of generalized inverse, it is necessary to study the novel TVMI solution models and find equivalence between existing solution models from PEs. But as far as we know, there is no relevant research so far. This paper focuses on the PE relevant models. By constructing different matrix error-monitoring functions from PEs, three corresponding time-variant matrix inverse solution models are derived. Through theoretical derivation and computer simulation, the feasibility and the effectiveness of the new models in TVMI problem solving are verified.The remainder of this paper is organized as follows. Section2 explains the PEs and some necessary equations. In Section 3, we derive a new ZNN multiple-multiplication model for matrix inversion, abbreviated as ZMMMI, from PEs. Thereafter, we employ the Euler forward finite difference (two-instant) formula, the four-instant Taylor–Zhang formula, the six-instant Zhang time discretization (ZTD) formula [31–33], the eight-instant ZTD formula, and the ten-instant ZTD formula to develop five discrete algorithms to compute TVMI. In Section 4, we derive a new PEs-based Getz–Marsden dynamic system (PGMDS) model for time-variant matrix inversion from PEs. Thereafter, we employ the five discrete formulas to present five discrete algorithms to compute TVMI. In Section 5, we present a DDD (direct derivative dynamics) neural network model for time-variant matrix inversion on the basis of PEs. Thereafter, we employ the five discrete formulas to develop five discrete algorithms to compute TVMI. In Section 6, we conduct numerical experiments to verify the convergence and precision of the three new models. Based on the experiment results, we compare and discuss the effectiveness of the three models. Section 7 gives the conclusions and future directions of research.Through the research of this paper, more TVMI solution models will be available to researchers. The main contributions of this paper can be outlined as follows.(1)
Three different models for solving the TVMI are proposed by defining different error-monitoring functions from PEs for the first time. This is the main motivation of this paper.(2)
We propose and provide a novel ZMMMI model from PEs, of which the stability and convergence are proved theoretically, thereby being rare complements to the problem solving of TVMI.(3)
The paper investigates and provides the theoretical analysis of the continuous-time models and finds that the classic GMDS model can also be derived from the PEs.(4)
For comparison, by exploiting the multiple-instant ZTD formulas, five discrete algorithms for each of the continuous-time ZMMMI, PGMDS, and DDD models are presented and constructed, respectively.(5)
The stability and convergence of the ZMMMI2i through ZMMMI10i algorithms and PGMDS2i through PGMDS10i algorithms are substantiated by numerical experiments of two benchmark examples.(6)
Numerical experiment results further indicate that compared with the effective ZMMMI algorithms and PGMDS algorithms, the DDD algorithms are less effective for TVMI.


## 2. Penrose Equations

For the purpose of laying a basis for direction, some necessary preliminaries of the time-variant matrix pseudoinverse are given. For a given real time-variant matrixAt∈ℝm×n, if one assumes that there exists an unknown real time-variant Xt∈ℝn×m which [27, 34](1a)AtXtAt=At,(1b)XtAtXt=Xt,(1c)AtXtT=AtXt,(1d)XtAtT=XtAt,where t denotes the time, superscript T indicates the matrix transpose operator, and Xt represents the time-variant pseudoinverse of At, which is often denoted by A+t. When rank (At) is equal to m and AtATt is nonsingular, the unique Xt can be obtained as(2)X+t≔ATtAtATt−1∈ℝn×m.Equation (2) holds true when At is full row rank. Similarly, when At is full column rank and AtATt is nonsingular, unique Xt can be obtained as(3)X+t≔ATtAt−1ATt∈ℝn×m.In equations (2) and (3), we use superscript −1 to indicate the matrix inversion operator. Also, X+t is called the right and left pseudoinverse of At in (2) and (3), respectively.


## 3. New ZNN Model and Algorithms

This section focuses on constructing another new ZF from PEs and then derives a new continuous-time solution model for Penrose pseudoinverse based on ZNN design formula [35].


### 3.1. Continuous-Time Model from PEs

The design formula of ZNN is shown below:(4)Z˙t=−γZt,where γ∈ℝ+ is the parameter whose physical meaning is the reciprocal of the product of the corresponding capacitance parameter and the resistance parameter. According to [28], (4) is a first-order differential equation whose norm of the general solution is inversely proportional to the exponential function of γ. Therefore, in order to make Zt converge to zero as rapid as possible, γ should be set as large as the physical hardware allows. Next, we discuss the ZNN model in the real situation. We choose (1a) to be the ZF as(5)Zt=AtXtAt−At,and when t⟶+∞, Zt⟶0 theoretically. Then, we take the derivative of both sides of (5) and get(6)Z˙t=A˙tXtAt+AtX˙tAt+AtXtA˙t−A˙t.Substituting (5) and (6) into (4), we obtain(7)Z˙t=−γAtXtAt−At=A˙tXtAt+AtX˙tAt+AtXtA˙t−A˙t.Reformulating (6), we have(8)AtX˙tAt=−A˙tXtAt−AtXtA˙t+A˙t−γAtXtAt−At.In order to derive the multiple-multiplication model, one assumes thatAt is square and of full rank. To get the explicit X˙t, we need to left multiply both sides of (8) by A−1t and right multiply both sides of (8) by A−1t and then get explicit X˙t as(9)X˙t=A−1t−A˙tXtAt−AtXtA˙t+A˙t−γAtXtAt−AtA−1t.The goal of the model design in this section is to use the known matrix to estimateA˙t, but A−1t appears on the right-hand side of (9), so it is necessary to use the substitution technique to replace A−1t with the appropriate known matrix. Based on Theorem 1 in [36], the state matrix Xt of (1a) globally converges to A−1t when the sampling gap ϵ⟶0 and t evolves large enough. Therefore, Xt would be a feasible substitute for A−1t. This substitution technique has been proven to be effective in [37–39]. Therefore, we have the explicit dynamics of the ZNN models as(10)X˙t=Xt−A˙tXtAt−AtXtA˙t+A˙t−γAtXtAt−AtXt.This novel model differs from any previous TVMI solution model. Note that the right-hand side of (10) involves multiple matrix-multiplication, so the model is termed as continuous-time ZMMMI model (10). To prepare for the following discussion, a lemma is given below [28].Lemma 1.
For a time-variant real matrixAt∈ℝn×n and its inverse A−1t, we have A˙−1t=−A−1tA˙tA−1t.Proof.
BecauseAtA−1t=I with I∈ℝn×n being the identity matrix, we have(11)dAtA−1tdt=dIdt=0.
Expanding the left-hand side of the above equation, we obtain(12)dAtdtA−1t+AtdA−1tdt=0.
Next, left multiplyA−1t on both sides of the above equation and then reformulate it. Thus, we obtain(13)dA−1tdt=−A−1tdAtdtA−1t.
For simplicity, the notationd⋅/dt is substituted by ⋅.. Then, we obtain(14)A˙−1t=−A−1tA˙tA−1t.
The proof is thus complete.For the continuous-time ZMMMI model (10), we provide the following proposition on its exponential convergence performance with a proper random initial state.Proposition 1.
For a smoothly time-variant real matrixAt∈ℝn×n of full rank, the state matrix Xt of ZMMMI model (10) starting from proper random initial state Xt exponentially converges to the time-variant theoretical inverse of At.Proof.
Letαt=Xt−A−1t represent the difference between the solution generated by ZMMMI model (10) and the theoretical inverse of At. Substituting Xt=αt+A−1t into (8), we obtain(15)Aα˙+A˙−1A=−A˙α+A−1A−Aα+A−1A˙+A˙−γAα+A−1A−A,where A, α, A˙, and α˙ denote the time-variant matrices At, αt, A˙t, and α˙t, respectively, for conciseness. According to Lemma 1 and A−1A=AA−1=I, (15) can be reformulated to(16)A˙αA+Aα˙A+AαA˙=−γAαA.
LetW=AαA. The above equation can be rewritten as(17)W˙=−γW,which is the compact matrix form of the following set of n×n first-order differential equations:(18)ω˙ijt=−γωijt,∀i,j∈1,2,…,n.
Evidently, a Lyapunov function candidatelij=ωij2/2≥0 can be defined for ijth subsystem (18), which is also positive definite, i.e., lij>0 for ωij≠0 and lij=0 for ωij=0. Therefore, we have its time derivative(19)dlijtdt=ωijω˙ij=−γωij2.
Evidently,l˙ij is negative definite, i.e., l˙ij<0 for ωij≠0 and l˙ij=0 for ωij=0. Besides, when ωij⟶∞, the Lyapunov function candidate lij=ωij2/2⟶∞, where symbol ⋅ represents the absolute value of a scalar value. According to the Lyapunov stability theory, ωijt (globally) converges to zero for any i,j∈1,2,…,n [40]. Therefore, for ZMMMI model (10), the state matrix Xt can converge to the time-variant theoretical inverse A−1t starting from proper random initial state X0, e.g., sufficiently close to A−10 [36, 41, 42]. Approximately and conditionally, we prove the exponential convergence performance of ZMMMI model (10).
We get the analytic solution of (18) in the following matrix form:(20)Wt=W0e−γt.
Then, we further obtain(21)WtF=W0Fe−γt,where symbol ⋅F denotes the Frobenius norm of a matrix and e represents Euler number. According to (21), as t⟶∞, WtF⟶0 with rate γ>0 exponentially. The proof of exponential convergence of ZMMMI model (10) is thus completed.


### 3.2. Discrete Algorithms

The hardware implementation of the continuous-time ZMMMI model (10) needs discrete algorithm. In this section, we discuss five discrete algorithms for the continuous-time ZMMMI model (10).


#### 3.2.1. Euler Forward Formula-Based Algorithm

For simplicity, the following Euler forward difference method is referred [43, 44]:(22)X˙tk=1ϵXtk+1−Xtk+Oϵ,where ϵ>0 represents the sampling gap and k=0,1,2,… represents the iteration number. tk=kϵ denotes the sampling time at iteration number k. For simplicity, we denote Xk=Xtk, and then (22) is transformed as(23)X˙k=1ϵXk+1−Xk+Oϵ.With (23), we discretize the continuous-time ZMMMI model (10) as(24)1ϵXk+1−Xk=Xk−A˙kXkAk−AkXkA˙k+A˙k−γAkXkAk−AkXk+Oϵ,which is further formulated as(25)Xk+1≐Xk+ϵXk−A˙kXkAk−AkXkA˙k+A˙kXk−hXkAkXkAk−AkXk,where h=ϵγ>0 denotes sufficiently small step length. Because it only depends on two sampling points in the past time, algorithm (25) is named as ZMMMI2i algorithm. It is worth pointing out that ϵ>0 should be set sufficiently small to ensure the convergence of the discrete model of ZMMMI model (10).


#### 3.2.2. Taylor–Zhang Discretization Formula-Based Algorithm

According to [45], TZDF is presented as(26)u˙k=1ϵuk+1−32uk+uk−1−12uk−2+Oϵ2.With (10) and (26), the four-instant ZMMMI algorithm is named as ZMMMI4i algorithm, which is expressed as follows:(27)Xk+1≐ϵXk−A˙kXkAk−AkXkA˙k+A˙kXk−hXkAkXkAk−AkXk+32Xk−Xk−1+12Xk−2,where h>0 denotes the step size again as before.


#### 3.2.3. Six-Instant ZTD Formula-Based Algorithm

In [37], a six-instant ZTD formula is given as follows:(28)u˙k=1ϵ12uk+1−548uk−14uk−1−18uk−2−112uk−3+116uk−4+Oϵ3.With (10) and (28), the six-instant ZMMMI algorithm is named as ZMMMI6i algorithm, which is given as follows:(29)Xk+1≐2ϵXk−A˙kXkAk−AkXkA˙k+A˙kXk−2hXkAkXkAk−AkXk+524Xk+12Xk−1+14Xk−2+16Xk−3−18Xk−4,where h>0 denotes the step size again as before.


#### 3.2.4. Eight-Instant ZTD Formula-Based Algorithm

In [37], an eight-instant ZTD formula is given as follows:(30)u˙k=1111ϵ50uk+1−5110uk−20uk−1−30uk−2−10uk−3+354uk−4+445uk−5−5uk−6+Oϵ4.With (10) and (30), the eight-instant ZMMMI algorithm is named as ZMMMI8i algorithm, which is given as follows:(31)Xk+1≐ϵXk−A˙kXkAk−AkXkA˙k+A˙kXk−hXkAkXkAk−AkXk+111505110Xk+20Xk−1+30Xk−2+10Xk−3−354Xk−4−445Xk−5+5Xk−6.


#### 3.2.5. Ten-Instant ZTD Formula-Based Algorithm

In [25], one ten-instant formula is given as follows:(32)u˙k=1ϵ∑i=09ciuk+1−i+Oϵ5,where c0=100/237, c1=947/66360, c2=−20/237, c3=−100/237, c4=−10/79, c5=40/237, c6=88/1185, c7=10/237, c8=−230/1659, and c9=95/1896 are the coefficients presented in [25]. With (10) and (32), the ten-instant ZMMMI algorithm is named as ZMMMI10i algorithm, which is thus obtained:(33)Xk+1≐1c0ϵXk−A˙kXkAk−AkXkA˙k+A˙kXk−hXkAkXkAk−AkXk−∑i=19ciXk+1−i,where h=ϵγ. It should be noticed that when h is a given constant, the truncation error is consistent with Oϵ6. When h is not fixed and γ is a constant, the error is upgraded to Oϵ5 as a result.


### 3.3. Steady-State Residual Errors of Discrete Algorithms

Next, a theorem is given to show that the steady-state residual errors oflimk⟶+∞supAk+1Xk+1Ak+1−Ak+1F are equivalent to the precision of the corresponding discrete algorithm in this section.Theorem 1.
Consider a smoothly time-variant real matrixAt∈ℝn×n of full rank. With sufficiently small sampling gap ϵ∈0,1, the maximal steady-state residual error limk⟶+∞supAk+1Xk+1Ak+1−Ak+1F of discrete algorithm (25), (27), (29), (31), or (33) is Oϵp, where p=2,3,4,5or6.Proof.
According to Theorems1 and 2 in [46], we know the fact that (25) is 0-stable, consistent, and convergent. As a result, it converges with the order of its truncation error. Assume that Bk+1 is the exact solution of Ak+1Bk+1=I. According to (25), (27), (29), (31), and (33), we have Xk+1=Bk+1+Oϵp, with ϵ∈0,1, and further have(34)=limk⟶+∞supAk+1Xk+1Ak+1−Ak+1F=limk⟶+∞supAk+1Bk+1+OϵpAk+1−Ak+1F=limk⟶+∞supI+Ak+1OϵpAk+1−Ak+1F=limk⟶+∞supAk+1+Ak+12Oϵp−Ak+1F=limk⟶+∞supAk+12OϵpF=Oϵp.
The proof is thus completed.


## 3.1. Continuous-Time Model from PEs

The design formula of ZNN is shown below:(4)Z˙t=−γZt,where γ∈ℝ+ is the parameter whose physical meaning is the reciprocal of the product of the corresponding capacitance parameter and the resistance parameter. According to [28], (4) is a first-order differential equation whose norm of the general solution is inversely proportional to the exponential function of γ. Therefore, in order to make Zt converge to zero as rapid as possible, γ should be set as large as the physical hardware allows. Next, we discuss the ZNN model in the real situation. We choose (1a) to be the ZF as(5)Zt=AtXtAt−At,and when t⟶+∞, Zt⟶0 theoretically. Then, we take the derivative of both sides of (5) and get(6)Z˙t=A˙tXtAt+AtX˙tAt+AtXtA˙t−A˙t.Substituting (5) and (6) into (4), we obtain(7)Z˙t=−γAtXtAt−At=A˙tXtAt+AtX˙tAt+AtXtA˙t−A˙t.Reformulating (6), we have(8)AtX˙tAt=−A˙tXtAt−AtXtA˙t+A˙t−γAtXtAt−At.In order to derive the multiple-multiplication model, one assumes thatAt is square and of full rank. To get the explicit X˙t, we need to left multiply both sides of (8) by A−1t and right multiply both sides of (8) by A−1t and then get explicit X˙t as(9)X˙t=A−1t−A˙tXtAt−AtXtA˙t+A˙t−γAtXtAt−AtA−1t.The goal of the model design in this section is to use the known matrix to estimateA˙t, but A−1t appears on the right-hand side of (9), so it is necessary to use the substitution technique to replace A−1t with the appropriate known matrix. Based on Theorem 1 in [36], the state matrix Xt of (1a) globally converges to A−1t when the sampling gap ϵ⟶0 and t evolves large enough. Therefore, Xt would be a feasible substitute for A−1t. This substitution technique has been proven to be effective in [37–39]. Therefore, we have the explicit dynamics of the ZNN models as(10)X˙t=Xt−A˙tXtAt−AtXtA˙t+A˙t−γAtXtAt−AtXt.This novel model differs from any previous TVMI solution model. Note that the right-hand side of (10) involves multiple matrix-multiplication, so the model is termed as continuous-time ZMMMI model (10). To prepare for the following discussion, a lemma is given below [28].Lemma 1.
For a time-variant real matrixAt∈ℝn×n and its inverse A−1t, we have A˙−1t=−A−1tA˙tA−1t.Proof.
BecauseAtA−1t=I with I∈ℝn×n being the identity matrix, we have(11)dAtA−1tdt=dIdt=0.
Expanding the left-hand side of the above equation, we obtain(12)dAtdtA−1t+AtdA−1tdt=0.
Next, left multiplyA−1t on both sides of the above equation and then reformulate it. Thus, we obtain(13)dA−1tdt=−A−1tdAtdtA−1t.
For simplicity, the notationd⋅/dt is substituted by ⋅.. Then, we obtain(14)A˙−1t=−A−1tA˙tA−1t.
The proof is thus complete.For the continuous-time ZMMMI model (10), we provide the following proposition on its exponential convergence performance with a proper random initial state.Proposition 1.
For a smoothly time-variant real matrixAt∈ℝn×n of full rank, the state matrix Xt of ZMMMI model (10) starting from proper random initial state Xt exponentially converges to the time-variant theoretical inverse of At.Proof.
Letαt=Xt−A−1t represent the difference between the solution generated by ZMMMI model (10) and the theoretical inverse of At. Substituting Xt=αt+A−1t into (8), we obtain(15)Aα˙+A˙−1A=−A˙α+A−1A−Aα+A−1A˙+A˙−γAα+A−1A−A,where A, α, A˙, and α˙ denote the time-variant matrices At, αt, A˙t, and α˙t, respectively, for conciseness. According to Lemma 1 and A−1A=AA−1=I, (15) can be reformulated to(16)A˙αA+Aα˙A+AαA˙=−γAαA.
LetW=AαA. The above equation can be rewritten as(17)W˙=−γW,which is the compact matrix form of the following set of n×n first-order differential equations:(18)ω˙ijt=−γωijt,∀i,j∈1,2,…,n.
Evidently, a Lyapunov function candidatelij=ωij2/2≥0 can be defined for ijth subsystem (18), which is also positive definite, i.e., lij>0 for ωij≠0 and lij=0 for ωij=0. Therefore, we have its time derivative(19)dlijtdt=ωijω˙ij=−γωij2.
Evidently,l˙ij is negative definite, i.e., l˙ij<0 for ωij≠0 and l˙ij=0 for ωij=0. Besides, when ωij⟶∞, the Lyapunov function candidate lij=ωij2/2⟶∞, where symbol ⋅ represents the absolute value of a scalar value. According to the Lyapunov stability theory, ωijt (globally) converges to zero for any i,j∈1,2,…,n [40]. Therefore, for ZMMMI model (10), the state matrix Xt can converge to the time-variant theoretical inverse A−1t starting from proper random initial state X0, e.g., sufficiently close to A−10 [36, 41, 42]. Approximately and conditionally, we prove the exponential convergence performance of ZMMMI model (10).
We get the analytic solution of (18) in the following matrix form:(20)Wt=W0e−γt.
Then, we further obtain(21)WtF=W0Fe−γt,where symbol ⋅F denotes the Frobenius norm of a matrix and e represents Euler number. According to (21), as t⟶∞, WtF⟶0 with rate γ>0 exponentially. The proof of exponential convergence of ZMMMI model (10) is thus completed.


## 3.2. Discrete Algorithms

The hardware implementation of the continuous-time ZMMMI model (10) needs discrete algorithm. In this section, we discuss five discrete algorithms for the continuous-time ZMMMI model (10).


### 3.2.1. Euler Forward Formula-Based Algorithm

For simplicity, the following Euler forward difference method is referred [43, 44]:(22)X˙tk=1ϵXtk+1−Xtk+Oϵ,where ϵ>0 represents the sampling gap and k=0,1,2,… represents the iteration number. tk=kϵ denotes the sampling time at iteration number k. For simplicity, we denote Xk=Xtk, and then (22) is transformed as(23)X˙k=1ϵXk+1−Xk+Oϵ.With (23), we discretize the continuous-time ZMMMI model (10) as(24)1ϵXk+1−Xk=Xk−A˙kXkAk−AkXkA˙k+A˙k−γAkXkAk−AkXk+Oϵ,which is further formulated as(25)Xk+1≐Xk+ϵXk−A˙kXkAk−AkXkA˙k+A˙kXk−hXkAkXkAk−AkXk,where h=ϵγ>0 denotes sufficiently small step length. Because it only depends on two sampling points in the past time, algorithm (25) is named as ZMMMI2i algorithm. It is worth pointing out that ϵ>0 should be set sufficiently small to ensure the convergence of the discrete model of ZMMMI model (10).


### 3.2.2. Taylor–Zhang Discretization Formula-Based Algorithm

According to [45], TZDF is presented as(26)u˙k=1ϵuk+1−32uk+uk−1−12uk−2+Oϵ2.With (10) and (26), the four-instant ZMMMI algorithm is named as ZMMMI4i algorithm, which is expressed as follows:(27)Xk+1≐ϵXk−A˙kXkAk−AkXkA˙k+A˙kXk−hXkAkXkAk−AkXk+32Xk−Xk−1+12Xk−2,where h>0 denotes the step size again as before.


### 3.2.3. Six-Instant ZTD Formula-Based Algorithm

In [37], a six-instant ZTD formula is given as follows:(28)u˙k=1ϵ12uk+1−548uk−14uk−1−18uk−2−112uk−3+116uk−4+Oϵ3.With (10) and (28), the six-instant ZMMMI algorithm is named as ZMMMI6i algorithm, which is given as follows:(29)Xk+1≐2ϵXk−A˙kXkAk−AkXkA˙k+A˙kXk−2hXkAkXkAk−AkXk+524Xk+12Xk−1+14Xk−2+16Xk−3−18Xk−4,where h>0 denotes the step size again as before.


### 3.2.4. Eight-Instant ZTD Formula-Based Algorithm

In [37], an eight-instant ZTD formula is given as follows:(30)u˙k=1111ϵ50uk+1−5110uk−20uk−1−30uk−2−10uk−3+354uk−4+445uk−5−5uk−6+Oϵ4.With (10) and (30), the eight-instant ZMMMI algorithm is named as ZMMMI8i algorithm, which is given as follows:(31)Xk+1≐ϵXk−A˙kXkAk−AkXkA˙k+A˙kXk−hXkAkXkAk−AkXk+111505110Xk+20Xk−1+30Xk−2+10Xk−3−354Xk−4−445Xk−5+5Xk−6.


### 3.2.5. Ten-Instant ZTD Formula-Based Algorithm

In [25], one ten-instant formula is given as follows:(32)u˙k=1ϵ∑i=09ciuk+1−i+Oϵ5,where c0=100/237, c1=947/66360, c2=−20/237, c3=−100/237, c4=−10/79, c5=40/237, c6=88/1185, c7=10/237, c8=−230/1659, and c9=95/1896 are the coefficients presented in [25]. With (10) and (32), the ten-instant ZMMMI algorithm is named as ZMMMI10i algorithm, which is thus obtained:(33)Xk+1≐1c0ϵXk−A˙kXkAk−AkXkA˙k+A˙kXk−hXkAkXkAk−AkXk−∑i=19ciXk+1−i,where h=ϵγ. It should be noticed that when h is a given constant, the truncation error is consistent with Oϵ6. When h is not fixed and γ is a constant, the error is upgraded to Oϵ5 as a result.


## 3.2.1. Euler Forward Formula-Based Algorithm

For simplicity, the following Euler forward difference method is referred [43, 44]:(22)X˙tk=1ϵXtk+1−Xtk+Oϵ,where ϵ>0 represents the sampling gap and k=0,1,2,… represents the iteration number. tk=kϵ denotes the sampling time at iteration number k. For simplicity, we denote Xk=Xtk, and then (22) is transformed as(23)X˙k=1ϵXk+1−Xk+Oϵ.With (23), we discretize the continuous-time ZMMMI model (10) as(24)1ϵXk+1−Xk=Xk−A˙kXkAk−AkXkA˙k+A˙k−γAkXkAk−AkXk+Oϵ,which is further formulated as(25)Xk+1≐Xk+ϵXk−A˙kXkAk−AkXkA˙k+A˙kXk−hXkAkXkAk−AkXk,where h=ϵγ>0 denotes sufficiently small step length. Because it only depends on two sampling points in the past time, algorithm (25) is named as ZMMMI2i algorithm. It is worth pointing out that ϵ>0 should be set sufficiently small to ensure the convergence of the discrete model of ZMMMI model (10).


## 3.2.2. Taylor–Zhang Discretization Formula-Based Algorithm

According to [45], TZDF is presented as(26)u˙k=1ϵuk+1−32uk+uk−1−12uk−2+Oϵ2.With (10) and (26), the four-instant ZMMMI algorithm is named as ZMMMI4i algorithm, which is expressed as follows:(27)Xk+1≐ϵXk−A˙kXkAk−AkXkA˙k+A˙kXk−hXkAkXkAk−AkXk+32Xk−Xk−1+12Xk−2,where h>0 denotes the step size again as before.


## 3.2.3. Six-Instant ZTD Formula-Based Algorithm

In [37], a six-instant ZTD formula is given as follows:(28)u˙k=1ϵ12uk+1−548uk−14uk−1−18uk−2−112uk−3+116uk−4+Oϵ3.With (10) and (28), the six-instant ZMMMI algorithm is named as ZMMMI6i algorithm, which is given as follows:(29)Xk+1≐2ϵXk−A˙kXkAk−AkXkA˙k+A˙kXk−2hXkAkXkAk−AkXk+524Xk+12Xk−1+14Xk−2+16Xk−3−18Xk−4,where h>0 denotes the step size again as before.


## 3.2.4. Eight-Instant ZTD Formula-Based Algorithm

In [37], an eight-instant ZTD formula is given as follows:(30)u˙k=1111ϵ50uk+1−5110uk−20uk−1−30uk−2−10uk−3+354uk−4+445uk−5−5uk−6+Oϵ4.With (10) and (30), the eight-instant ZMMMI algorithm is named as ZMMMI8i algorithm, which is given as follows:(31)Xk+1≐ϵXk−A˙kXkAk−AkXkA˙k+A˙kXk−hXkAkXkAk−AkXk+111505110Xk+20Xk−1+30Xk−2+10Xk−3−354Xk−4−445Xk−5+5Xk−6.


## 3.2.5. Ten-Instant ZTD Formula-Based Algorithm

In [25], one ten-instant formula is given as follows:(32)u˙k=1ϵ∑i=09ciuk+1−i+Oϵ5,where c0=100/237, c1=947/66360, c2=−20/237, c3=−100/237, c4=−10/79, c5=40/237, c6=88/1185, c7=10/237, c8=−230/1659, and c9=95/1896 are the coefficients presented in [25]. With (10) and (32), the ten-instant ZMMMI algorithm is named as ZMMMI10i algorithm, which is thus obtained:(33)Xk+1≐1c0ϵXk−A˙kXkAk−AkXkA˙k+A˙kXk−hXkAkXkAk−AkXk−∑i=19ciXk+1−i,where h=ϵγ. It should be noticed that when h is a given constant, the truncation error is consistent with Oϵ6. When h is not fixed and γ is a constant, the error is upgraded to Oϵ5 as a result.


## 3.3. Steady-State Residual Errors of Discrete Algorithms

Next, a theorem is given to show that the steady-state residual errors oflimk⟶+∞supAk+1Xk+1Ak+1−Ak+1F are equivalent to the precision of the corresponding discrete algorithm in this section.Theorem 1.
Consider a smoothly time-variant real matrixAt∈ℝn×n of full rank. With sufficiently small sampling gap ϵ∈0,1, the maximal steady-state residual error limk⟶+∞supAk+1Xk+1Ak+1−Ak+1F of discrete algorithm (25), (27), (29), (31), or (33) is Oϵp, where p=2,3,4,5or6.Proof.
According to Theorems1 and 2 in [46], we know the fact that (25) is 0-stable, consistent, and convergent. As a result, it converges with the order of its truncation error. Assume that Bk+1 is the exact solution of Ak+1Bk+1=I. According to (25), (27), (29), (31), and (33), we have Xk+1=Bk+1+Oϵp, with ϵ∈0,1, and further have(34)=limk⟶+∞supAk+1Xk+1Ak+1−Ak+1F=limk⟶+∞supAk+1Bk+1+OϵpAk+1−Ak+1F=limk⟶+∞supI+Ak+1OϵpAk+1−Ak+1F=limk⟶+∞supAk+1+Ak+12Oϵp−Ak+1F=limk⟶+∞supAk+12OϵpF=Oϵp.
The proof is thus completed.


## 4. GMDS Model and Algorithms

This section focuses on constructing new ZF from PEs and then derives continuous-time PGMDS solution model for matrix inverse based on ZNN design formula. Furthermore, we develop five discrete algorithms for continuous-time PGMDS model by exploiting the multiple-instant ZTD formulas.


### 4.1. Continuous-Time Model from PEs

In this section, we choose (1b) to be the ZF as(35)Zt=XtAtXt−Xt,and when t⟶+∞, Zt⟶0 theoretically. Then, we take the derivative of both sides of (1b) and get(36)Z˙t=X˙tAtXt+XtA˙tXt+XtAtX˙t−X˙t.Substituting (35) and (36) into (4), we obtain(37)X˙tAtXt+XtA˙tXt+XtAtX˙t−X˙t=−γXtAtXt−Xt.According to the fact thatAtXt=XtAt=I, we reformulate (36) and then obtain(38)X˙t=−XtA˙tXt−γXtAtXt−Xt.Model (38), which is derived from PE (1b), is exactly the Getz–Marsden dynamic system mentioned in [36]. Therefore, the model is named as continuous-time PGMDS model (38).


### 4.2. Discrete Algorithms

In this part, we discretize the continuous-time PGMDS model (38) with the following five discrete formulas as those in Section 3.


#### 4.2.1. Euler Forward Formula-Based Algorithm

First, (38) is discretized as(39)X˙k=−XkA˙kXk−γXkAkXk−Xk+Oϵ.With (39) and (23), the two-instant PGMDS algorithm is named as PGMDS2i algorithm, which is expressed as follows:(40)Xk+1≐Xk−ϵXkA˙kXk−hXkAkXk−Xk,where h=ϵγ>0 denotes the step length as before.


#### 4.2.2. Taylor–Zhang Discretization Formula-Based Algorithm

With (39) and (26), the four-instant PGMDS algorithm is named as PGMDS4i algorithm, which is expressed as follows:(41)Xk+1≐−ϵXkA˙kXk−hXkAkXk−Xk+32Xk−Xk−1+12Xk−2,where h>0 denotes the step size again as before.


#### 4.2.3. Six-Instant ZTD Formula-Based Algorithm

With (39) and (28), the six-instant PGMDS algorithm is named as PGMDS6i algorithm, which is given as follows:(42)Xk+1≐−ϵ2XkA˙kXk−2hXkAkXk−Xk+524Xk+12Xk−1+14Xk−2+16Xk−3−18Xk−4,where h>0 denotes the step size again as before.


#### 4.2.4. Eight-Instant ZTD Formula-Based Algorithm

With (39) and (30), the eight-instant PGMDS algorithm is named as PGMDS8i algorithm, which is given as follows:(43)Xk+1≐−ϵXkA˙kXk−hXkAkXk−Xk+111505110Xk+20Xk−1+30Xk−2+10Xk−3−354Xk−4−445Xk−5+5Xk−6,where h=ϵγ>0 denotes the step length as before.


#### 4.2.5. Ten-Instant ZTD Formula-Based Algorithm

With (39) and (32), the ten-instant PGMDS algorithm is named as PGMDS10i algorithm, which is thus obtained:(44)Xk+1≐1c0−ϵXkA˙kXk−hXkAkXk−Xk−∑i=19ciXk+1−i,where h=ϵγ>0 denotes the step length as before.


### 4.3. Steady-State Residual Errors of Discrete Algorithms

Next, a theorem is given to show that the steady-state residual oflimk⟶+∞supXk+1Ak+1Xk+1−Xk+1F is equivalent to the precision of the corresponding discrete algorithm in this section.Theorem 2.
Consider a smoothly time-variant real matrixAt∈ℝn×n of full rank. With sufficiently small sampling gap ϵ∈0,1, the maximal steady-state residual error limk⟶+∞supXk+1Ak+1Xk+1−Xk+1F of discrete algorithm (40)–(44) is Oϵp, where p=2,3,4,5,and6, respectively.Proof.
With the same conditions as Theorem1, we know the fact that (39) is 0-stable, consistent, and convergent [46]. As a result, it converges with the order of its truncation error. Assume that Bk+1 satisfies Ak+1Bk+1=I and Xk+1=Bk+1+Oϵp with ϵ∈0,1 as before. We further have(45)=limk⟶+∞supXk+1Ak+1Xk+1−Xk+1F=limk⟶+∞supBk+1+OϵpAk+1Xk+1−Xk+1F=limk⟶+∞supBk+1Ak+1+Ak+1OϵpXk+1−Xk+1F=limk⟶+∞supI+Ak+1OϵpXk+1−Xk+1F=limk⟶+∞supXk+1+Ak+1Xk+1Oϵp−Xk+1F=limk⟶+∞supAk+1Bk+1+OϵpOϵpF=limk⟶+∞supOϵp+Ak+1Oϵ2pF=Oϵp.
The proof is thus completed.


## 4.1. Continuous-Time Model from PEs

In this section, we choose (1b) to be the ZF as(35)Zt=XtAtXt−Xt,and when t⟶+∞, Zt⟶0 theoretically. Then, we take the derivative of both sides of (1b) and get(36)Z˙t=X˙tAtXt+XtA˙tXt+XtAtX˙t−X˙t.Substituting (35) and (36) into (4), we obtain(37)X˙tAtXt+XtA˙tXt+XtAtX˙t−X˙t=−γXtAtXt−Xt.According to the fact thatAtXt=XtAt=I, we reformulate (36) and then obtain(38)X˙t=−XtA˙tXt−γXtAtXt−Xt.Model (38), which is derived from PE (1b), is exactly the Getz–Marsden dynamic system mentioned in [36]. Therefore, the model is named as continuous-time PGMDS model (38).


## 4.2. Discrete Algorithms

In this part, we discretize the continuous-time PGMDS model (38) with the following five discrete formulas as those in Section 3.


### 4.2.1. Euler Forward Formula-Based Algorithm

First, (38) is discretized as(39)X˙k=−XkA˙kXk−γXkAkXk−Xk+Oϵ.With (39) and (23), the two-instant PGMDS algorithm is named as PGMDS2i algorithm, which is expressed as follows:(40)Xk+1≐Xk−ϵXkA˙kXk−hXkAkXk−Xk,where h=ϵγ>0 denotes the step length as before.


### 4.2.2. Taylor–Zhang Discretization Formula-Based Algorithm

With (39) and (26), the four-instant PGMDS algorithm is named as PGMDS4i algorithm, which is expressed as follows:(41)Xk+1≐−ϵXkA˙kXk−hXkAkXk−Xk+32Xk−Xk−1+12Xk−2,where h>0 denotes the step size again as before.


### 4.2.3. Six-Instant ZTD Formula-Based Algorithm

With (39) and (28), the six-instant PGMDS algorithm is named as PGMDS6i algorithm, which is given as follows:(42)Xk+1≐−ϵ2XkA˙kXk−2hXkAkXk−Xk+524Xk+12Xk−1+14Xk−2+16Xk−3−18Xk−4,where h>0 denotes the step size again as before.


### 4.2.4. Eight-Instant ZTD Formula-Based Algorithm

With (39) and (30), the eight-instant PGMDS algorithm is named as PGMDS8i algorithm, which is given as follows:(43)Xk+1≐−ϵXkA˙kXk−hXkAkXk−Xk+111505110Xk+20Xk−1+30Xk−2+10Xk−3−354Xk−4−445Xk−5+5Xk−6,where h=ϵγ>0 denotes the step length as before.


### 4.2.5. Ten-Instant ZTD Formula-Based Algorithm

With (39) and (32), the ten-instant PGMDS algorithm is named as PGMDS10i algorithm, which is thus obtained:(44)Xk+1≐1c0−ϵXkA˙kXk−hXkAkXk−Xk−∑i=19ciXk+1−i,where h=ϵγ>0 denotes the step length as before.


## 4.2.1. Euler Forward Formula-Based Algorithm

First, (38) is discretized as(39)X˙k=−XkA˙kXk−γXkAkXk−Xk+Oϵ.With (39) and (23), the two-instant PGMDS algorithm is named as PGMDS2i algorithm, which is expressed as follows:(40)Xk+1≐Xk−ϵXkA˙kXk−hXkAkXk−Xk,where h=ϵγ>0 denotes the step length as before.


## 4.2.2. Taylor–Zhang Discretization Formula-Based Algorithm

With (39) and (26), the four-instant PGMDS algorithm is named as PGMDS4i algorithm, which is expressed as follows:(41)Xk+1≐−ϵXkA˙kXk−hXkAkXk−Xk+32Xk−Xk−1+12Xk−2,where h>0 denotes the step size again as before.


## 4.2.3. Six-Instant ZTD Formula-Based Algorithm

With (39) and (28), the six-instant PGMDS algorithm is named as PGMDS6i algorithm, which is given as follows:(42)Xk+1≐−ϵ2XkA˙kXk−2hXkAkXk−Xk+524Xk+12Xk−1+14Xk−2+16Xk−3−18Xk−4,where h>0 denotes the step size again as before.


## 4.2.4. Eight-Instant ZTD Formula-Based Algorithm

With (39) and (30), the eight-instant PGMDS algorithm is named as PGMDS8i algorithm, which is given as follows:(43)Xk+1≐−ϵXkA˙kXk−hXkAkXk−Xk+111505110Xk+20Xk−1+30Xk−2+10Xk−3−354Xk−4−445Xk−5+5Xk−6,where h=ϵγ>0 denotes the step length as before.


## 4.2.5. Ten-Instant ZTD Formula-Based Algorithm

With (39) and (32), the ten-instant PGMDS algorithm is named as PGMDS10i algorithm, which is thus obtained:(44)Xk+1≐1c0−ϵXkA˙kXk−hXkAkXk−Xk−∑i=19ciXk+1−i,where h=ϵγ>0 denotes the step length as before.


## 4.3. Steady-State Residual Errors of Discrete Algorithms

Next, a theorem is given to show that the steady-state residual oflimk⟶+∞supXk+1Ak+1Xk+1−Xk+1F is equivalent to the precision of the corresponding discrete algorithm in this section.Theorem 2.
Consider a smoothly time-variant real matrixAt∈ℝn×n of full rank. With sufficiently small sampling gap ϵ∈0,1, the maximal steady-state residual error limk⟶+∞supXk+1Ak+1Xk+1−Xk+1F of discrete algorithm (40)–(44) is Oϵp, where p=2,3,4,5,and6, respectively.Proof.
With the same conditions as Theorem1, we know the fact that (39) is 0-stable, consistent, and convergent [46]. As a result, it converges with the order of its truncation error. Assume that Bk+1 satisfies Ak+1Bk+1=I and Xk+1=Bk+1+Oϵp with ϵ∈0,1 as before. We further have(45)=limk⟶+∞supXk+1Ak+1Xk+1−Xk+1F=limk⟶+∞supBk+1+OϵpAk+1Xk+1−Xk+1F=limk⟶+∞supBk+1Ak+1+Ak+1OϵpXk+1−Xk+1F=limk⟶+∞supI+Ak+1OϵpXk+1−Xk+1F=limk⟶+∞supXk+1+Ak+1Xk+1Oϵp−Xk+1F=limk⟶+∞supAk+1Bk+1+OϵpOϵpF=limk⟶+∞supOϵp+Ak+1Oϵ2pF=Oϵp.
The proof is thus completed.


## 5. DDD Model and Algorithms

This section focuses on constructing direct derivative dynamics from one appropriate Penrose equation and then derives direct derivative solution model for matrix inversion.


### 5.1. Continuous-Time Model from PEs

We take the derivative on both sides of (1a) and then obtain(46)A˙tXtAt+AtX˙tAt+AtXtA˙t=A˙t.To get the explicit expression ofX˙t, we reformulate above equation and obtain(47)X˙t=−A−1tA˙tXt−XtA˙tA−1t+A−1tA˙tA−1t.According to the assumption of the matrix inverse problem, we can substituteA−1t with Xt in (47) and thus obtain(48)X˙t=−XtA˙tXt−XtA˙tXt+XtA˙tXt.The right-hand side of above equation can be simplified as(49)X˙t=−XtA˙tXt.Finally, we present the continuous-time direct derivative dynamics model (49), which is termed as continuous-time DDD model (49) for short.


### 5.2. Discrete Algorithms

In this part, we discretize the continuous-time DDD model (49) with the following five discrete formulas as those in Section 3.


#### 5.2.1. Euler Forward Formula-Based Algorithm

First, (49) is discretized as(50)X˙k=−XkAk.Xk+Oϵ.With (50) and (23), the two-instant DDD algorithm is named as DDD2i algorithm, which is expressed as follows:(51)Xk+1≐Xk−ϵXkA˙Xk.


#### 5.2.2. Taylor–Zhang Discretization Formula-Based Algorithm

With (50) and (26), the four-instant DDD algorithm is named as DDD4i algorithm, which is expressed as follows:(52)Xk+1≐−ϵXkAk.Xk+32Xk−Xk−1+12Xk−2.


#### 5.2.3. Six-Instant ZTD Formula-Based Algorithm

With (50) and (28), the six-instant DDD algorithm is named as DDD6i algorithm, which is given as follows:(53)Xk+1≐−ϵXkAk.Xk+524Xk+12Xk−1+14Xk−2+16Xk−3−18Xk−4.


#### 5.2.4. Eight-Instant ZTD Formula-Based Algorithm

With (50) and (30), the eight-instant DDD algorithm is named as DDD8i algorithm, which is given as follows:(54)Xk+1≐−ϵXkAk.Xk+111505110Xk+20Xk−1+30Xk−2+10Xk−3−354Xk−4−445Xk−5+5Xk−6.


#### 5.2.5. Ten-Instant ZTD Formula-Based Algorithm

With (50) and (32), the ten-instant DDD algorithm is named as DDD10i algorithm, which is thus obtained:(55)Xk+1≐−ϵXkAk.Xk−∑i=19ciXk+1−i.


## 5.1. Continuous-Time Model from PEs

We take the derivative on both sides of (1a) and then obtain(46)A˙tXtAt+AtX˙tAt+AtXtA˙t=A˙t.To get the explicit expression ofX˙t, we reformulate above equation and obtain(47)X˙t=−A−1tA˙tXt−XtA˙tA−1t+A−1tA˙tA−1t.According to the assumption of the matrix inverse problem, we can substituteA−1t with Xt in (47) and thus obtain(48)X˙t=−XtA˙tXt−XtA˙tXt+XtA˙tXt.The right-hand side of above equation can be simplified as(49)X˙t=−XtA˙tXt.Finally, we present the continuous-time direct derivative dynamics model (49), which is termed as continuous-time DDD model (49) for short.


## 5.2. Discrete Algorithms

In this part, we discretize the continuous-time DDD model (49) with the following five discrete formulas as those in Section 3.


### 5.2.1. Euler Forward Formula-Based Algorithm

First, (49) is discretized as(50)X˙k=−XkAk.Xk+Oϵ.With (50) and (23), the two-instant DDD algorithm is named as DDD2i algorithm, which is expressed as follows:(51)Xk+1≐Xk−ϵXkA˙Xk.


### 5.2.2. Taylor–Zhang Discretization Formula-Based Algorithm

With (50) and (26), the four-instant DDD algorithm is named as DDD4i algorithm, which is expressed as follows:(52)Xk+1≐−ϵXkAk.Xk+32Xk−Xk−1+12Xk−2.


### 5.2.3. Six-Instant ZTD Formula-Based Algorithm

With (50) and (28), the six-instant DDD algorithm is named as DDD6i algorithm, which is given as follows:(53)Xk+1≐−ϵXkAk.Xk+524Xk+12Xk−1+14Xk−2+16Xk−3−18Xk−4.


### 5.2.4. Eight-Instant ZTD Formula-Based Algorithm

With (50) and (30), the eight-instant DDD algorithm is named as DDD8i algorithm, which is given as follows:(54)Xk+1≐−ϵXkAk.Xk+111505110Xk+20Xk−1+30Xk−2+10Xk−3−354Xk−4−445Xk−5+5Xk−6.


### 5.2.5. Ten-Instant ZTD Formula-Based Algorithm

With (50) and (32), the ten-instant DDD algorithm is named as DDD10i algorithm, which is thus obtained:(55)Xk+1≐−ϵXkAk.Xk−∑i=19ciXk+1−i.


## 5.2.1. Euler Forward Formula-Based Algorithm

First, (49) is discretized as(50)X˙k=−XkAk.Xk+Oϵ.With (50) and (23), the two-instant DDD algorithm is named as DDD2i algorithm, which is expressed as follows:(51)Xk+1≐Xk−ϵXkA˙Xk.


## 5.2.2. Taylor–Zhang Discretization Formula-Based Algorithm

With (50) and (26), the four-instant DDD algorithm is named as DDD4i algorithm, which is expressed as follows:(52)Xk+1≐−ϵXkAk.Xk+32Xk−Xk−1+12Xk−2.


## 5.2.3. Six-Instant ZTD Formula-Based Algorithm

With (50) and (28), the six-instant DDD algorithm is named as DDD6i algorithm, which is given as follows:(53)Xk+1≐−ϵXkAk.Xk+524Xk+12Xk−1+14Xk−2+16Xk−3−18Xk−4.


## 5.2.4. Eight-Instant ZTD Formula-Based Algorithm

With (50) and (30), the eight-instant DDD algorithm is named as DDD8i algorithm, which is given as follows:(54)Xk+1≐−ϵXkAk.Xk+111505110Xk+20Xk−1+30Xk−2+10Xk−3−354Xk−4−445Xk−5+5Xk−6.


## 5.2.5. Ten-Instant ZTD Formula-Based Algorithm

With (50) and (32), the ten-instant DDD algorithm is named as DDD10i algorithm, which is thus obtained:(55)Xk+1≐−ϵXkAk.Xk−∑i=19ciXk+1−i.


## 6. Computer Experiments and Results

In this section, computer experiments is carried out to verify the effectiveness of the presented three models, ZMMMI model (10), PGMDS model (38), and DDD model (49), on three time-variant matrix inversion examples.


### 6.1. Example 1

Let us consider the following discrete-time matrix inversion problem withXk+1 to be obtained during tk,tk+1, of which Ak is defined as(56)Atk=sintk−costkcostksintk∈ℝ2×2.The task duration (i.e., final time) is uniformly set astd=30 s. To verify the computational results, the theoretical inversion of matrix (56) can be obtained as(57)Ytk=A−1tk=sintkcostk−costksintk∈ℝ2×2,which is given for checking the effectiveness and correctness of the three models and corresponding algorithms. The initial values of the example are arbitrarily set as(58)Y0=1.51.92.35.8,for all algorithms.


#### 6.1.1. New Models of Zhang Neural Network

Firstly, we check the effectiveness and the correctness of ZMMMI model (10) with five discrete algorithms, i.e., ZMMMI2i algorithm (25), ZMMMI4i algorithm (27), ZMMMI6i algorithm (29), ZMMMI8i algorithm (31), and ZMMMI10i algorithm (33). The step size and task duration are uniformly set as h=0.005 and td=30 s. The corresponding results of computer experiments are shown in Figures 1 and 2. Figure 1(a) shows the residual errors Xk+1Ak+1−IF of the five ZMMMI algorithms, with ϵ=0.001 s and h=0.005. Figure 1(b) illustrates the residual errors Xk+1−Ak+1−1F of the same five ZMMMI algorithms as before. From the trajectories of all entries in Figures 2(a)–2(d), the solution of the model coincides with the theoretical solution perfectly. In addition, we see that the residual error trajectories synthesized by different ZMMMI algorithms quickly stabilized to the steady-state error level, after undergoing the initial hundreds of recursions.Figure 1
Residual error trajectories synthesized by five ZMMMI algorithms for example 1 (withϵ=0.001 s and γ=5.0). (a) Trajectories of Xk+1Ak+1−IF. (b) Trajectories of Xk+1−Ak+1−1F.
(a)(b)Figure 2
Solution trajectories synthesized byYk+1 (i.e., Ak+1−1) and by ZMMMI10i algorithm (33) for example 1 (with ϵ=0.001 s and γ=5.0). (a) Trajectories of y11 and x11. (b) Trajectories of y12 and x12. (c) Trajectories of y21 and x21. (d) Trajectories of y22 and x22.
(a)(b)(c)(d)


#### 6.1.2. PGMDS Model

First, we check the effectiveness and the correctness of PGMDS model (38) with five discrete algorithms, i.e., PGMDS2i algorithm (40), PGMDS4i algorithm (41), PGMDS6i algorithm (42), PGMDS8i algorithm (43), and PGMDS10i algorithm (44). The step size and task duration are uniformly set as h=0.005 and td=30 s. The corresponding results of computer experiments are shown in Figures 3 and 4. Figure 3(a) shows the residual errors Xk+1Ak+1−IF of the five PGMDS algorithms, with ϵ=0.001 s and h=0.005. Figure 3(b) illustrates the residual errors Xk+1−Ak+1−1F of the same five PGMDS algorithms as before. From the trajectories of all entries in Figures 4(a)–4(d), the solution of the model coincides with the theoretical solution perfectly. In addition, we see that the residual error trajectories synthesized by different PGMDS algorithms quickly stabilized to the theoretical error level, after undergoing the initial hundreds of recursions.Figure 3
Residual error trajectories synthesized by five PGMDS algorithms for example 1 (withϵ=0.001 s and γ=5.0). (a) Trajectories of Xk+1Ak+1−IF. (b) Trajectories of Xk+1−Ak+1−1F.
(a)(b)Figure 4
Solution trajectories synthesized byYk+1 (i.e., Ak+1−1) and by PGMDS10i algorithm (39) for example 1 (with ϵ=0.001 s and γ=5.0). (a) Trajectories of y11 and x11. (b) Trajectories of y12 and x12. (c) Trajectories of y21 and x21. (d) Trajectories of y22 and x22.
(a)(b)(c)(d)


#### 6.1.3. DDD Model

First, we check the effectiveness and correctness of DDD model (49) with five discrete algorithms, i.e., DDD2i algorithm (51), DDD4i algorithm (52), DDD6i algorithm (53), DDD8i algorithm (54), and DDD10i algorithm (55). The step size and task duration are uniformly set as h=0.0001 and td=30 s. The corresponding results of computer experiments are shown in Figures 5 and 6. Figure 5(a) shows the residual errors Xk+1Ak+1−IF of the five DDD algorithms, with ϵ=0.001 s and h=0.005. Figure 5(b) illustrates the residual errors Xk+1−Ak+1−1F of the same five DDD algorithms as before. From the trajectories of all entries in Figures 6(a) and 6(b), even for the model with the highest accuracy, DDD10i algorithm (55), its solution still cannot converge to the theoretical solution.Figure 5
Residual error trajectories synthesized by five DDD algorithms for example 1 (withϵ=0.001 s and γ=5.0). (a) Trajectories of Xk+1Ak+1−IF. (b) Trajectories of Xk+1−Ak+1−1F.
(a)(b)Figure 6
Solution trajectories synthesized byYk+1 (i.e., Ak+1−1) and by DDD10i algorithm (55) for example 1 (with ϵ=0.001 s and γ=5.0). (a) Trajectories of y11 and x11. (b) Trajectories of y12 and x12. (c) Trajectories of y21 and x21. (d) Trajectories of y22 and x22.
(a)(b)(c)(d)


### 6.2. Example 2

The second time-variant matrix is a3×3 real matrix which is shown as follows:(59)Atk=5+costkcostk−sintkcostk3+sintksintkcostk−costk2+sintk,where Atk∈ℝ3×3. To verify the computational results, we utilize the theoretical inverse of Atk, which is denoted as Ytk=yijtk. Because Ytk is too complicated, it is omitted here. The initial values of the example are arbitrarily set as(60)Y0=9.50.70.6−1.1121.3−1.60.29.0.The task duration is uniformly set astd=30 s.


#### 6.2.1. ZMMMI Model

The results of the numerical experiments for ZMMMI2i algorithm (25), ZMMMI4i algorithm (27), ZMMMI6i algorithm (29), ZMMMI8i algorithm (31), and ZMMMI10i algorithm (33) are shown in Figures 7 and 8. The experimental results are satisfactory as expected.Figure 7
Residual error trajectories synthesized by five ZMMMI algorithms for example 2 (withϵ=0.0001 s and γ=2.0). (a) Trajectories of Xk+1Ak+1−IF. (b) Trajectories of Xk+1−Ak+1−1F.
(a)(b)Figure 8
Solution trajectories synthesized byYk+1 (i.e., Ak+1−1) and by ZMMMI10i algorithm (33) for example 2 (with ϵ=0.0001 s and γ=2.0). (a) Trajectories of y11 and x11. (b) Trajectories of y12 and x12. (c) Trajectories of y13 and x13. (d) Trajectories of y21 and x21. (e) Trajectories of y22 and x22. (f) Trajectories of y23 and x23. (g) Trajectories of y31 and x31. (h) Trajectories of y32 and x32. (i) Trajectories of y33 and x33.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


#### 6.2.2. PGMDS Model

The results of the numerical experiments for PGMDS2i algorithm (40), PGMDS4i algorithm (41), PGMDS6i algorithm (42), PGMDS8i algorithm (43), and PGMDS10i algorithm (44) are illustrated in Figures 9 and 10. The experimental results are satisfactory as expected.Figure 9
Residual error trajectories synthesized by five PGMDS algorithms for example 2 (withϵ=0.001 s and γ=5.0). (a) Trajectories of Xk+1Ak+1−IF. (b) Trajectories of Xk+1−Ak+1−1F.
(a)(b)Figure 10
Solution trajectories synthesized byYk+1 (i.e., Ak+1−1) and by PGMDS10i algorithm (44) for example 2 (with ϵ=0.001 s and γ=5.0). (a) Trajectories of y11 and x11. (b) Trajectories of y12 and x12. (c) Trajectories of y13 and x13. (d) Trajectories of y21 and x21. (e) Trajectories of y22 and x22. (f) Trajectories of y23 and x23. (g) Trajectories of y31 and x31. (h) Trajectories of y32 and x32. (i) Trajectories of y33 and x33.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


#### 6.2.3. DDD Model

The results of the numerical experiments for DDD2i algorithm (51), DDD4i algorithm (52), DDD6i algorithm (53), DDD8i algorithm (54), and DDD10i algorithm (55) are shown in Figures 11 and 12. In terms of convergence and accuracy, the experimental results in this example are consistent with those in example 1.Figure 11
Residual error trajectories synthesized by five DDD algorithms for example 2 (withϵ=0.001 s and γ=1.0). (a) Trajectories of Xk+1Ak+1−IF. (b) Trajectories of Xk+1−Ak+1−1F.
(a)(b)Figure 12
Solution trajectories synthesized byYk+1 (i.e., Ak+1−1) and by DDD10i algorithm (55) for example 2 (with ϵ=0.001 s and γ=1.0). (a) Trajectories of y11 and x11. (b) Trajectories of y12 and x12. (c) Trajectories of y13 and x13. (d) Trajectories of y21 and x21. (e) Trajectories of y22 and x22. (f) Trajectories of y23 and x23. (g) Trajectories of y31 and x31. (h) Trajectories of y32 and x32. (i) Trajectories of y33 and x33.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


### 6.3. Example 3

This example compares solely the effective models, i.e., ZMMMI (10) and PGMDS (38), with a 4×4 real matrix, which is shown as follows:(61)Atk=6+ckckckckck5+sk−skckckck4+sk−skskskck2+sk,where Atk∈ℝ4×4 and sk and ck denote sintk and costk, respectively. To verify the computational results, we utilize the theoretical inverse of Atk, which is denoted as Ytk=yijtk. Because Ytk is too complicated, it is omitted here. The initial values of the example are arbitrarily set as(62)Y0=6111150111400012.The task duration is uniformly set astd=30 s.


#### 6.3.1. ZMMMI Model

The results of the numerical experiments for ZMMMI2i algorithm (25), ZMMMI4i algorithm (27), ZMMMI6i algorithm (29), ZMMMI8i algorithm (31), and ZMMMI10i algorithm (33) are shown in Figures 13(a) and 13(b). The experimental results are satisfactory as expected.Figure 13
Residual error trajectories synthesized by five ZMMMI algorithms for example 3 (withϵ=0.0001 s and γ=2.0). (a) Trajectories of Xk+1Ak+1−IF. (b) Trajectories of Xk+1−Ak+1−1F.
(a)(b)


#### 6.3.2. PGMDS Model

The results of the numerical experiments for PGMDS2i algorithm (40), PGMDS4i algorithm (41), PGMDS6i algorithm (42), PGMDS8i algorithm (43), and PGMDS10i algorithm (44) are illustrated in Figures 14(a) and 14(b). The experimental results are satisfactory as expected.Figure 14
Residual error trajectories synthesized by five PGMDS algorithms for example 3 (withϵ=0.002 s and γ=5.0). (a) Trajectories of Xk+1Ak+1−IF. (b) Trajectories of Xk+1−Ak+1−1F.
(a)(b)


#### 6.3.3. Comparison and Discussion

We mainly investigate the residual error of the models and the coincidence between the solution matrix with the ground-truth matrix inverse. From the discrete simulation results of two examples of each five discrete algorithms, it can be seen that the convergence of ZMMMI model (10) and PGMDS model (38) is good, which is completely consistent with the conclusion of Proposition 1 and Theorems 1 and 2. Therefore, it can be concluded that the two models and corresponding discrete algorithms are effective. However, in the case of DDD model (49), the residual errors and the coincidence of the solution matrix entries are not satisfactory, so we evidently summarize that the effectiveness of DDD model (49) and corresponding discrete algorithms is not enough.


#### 6.3.4. Remark

Note that the convergence of DDD models is shown in [12, 47], where the DDD models are utilized to solve the time-varying nonlinear optimization problems. There are three points worth further discussing. First, in [12], the initial value is set to be the theoretical value, and the computational process is not disturbed. Second, comparatively, in this study, the initial value is not close enough to the theoretical value. Third, the higher-order discrete algorithm (i.e., the ten-instant discrete algorithm) is more easily affected by the rounding error disturbance [44]. Thus, the experimental results of the DDD model in this study showing the divergence are actually complementing the previous research studies, in addition to the confirmation of [47] about divergence. Therefore, we summarize that the DDD model is generally less effective, with further in-depth investigation being also a future research direction.


## 6.1. Example 1

Let us consider the following discrete-time matrix inversion problem withXk+1 to be obtained during tk,tk+1, of which Ak is defined as(56)Atk=sintk−costkcostksintk∈ℝ2×2.The task duration (i.e., final time) is uniformly set astd=30 s. To verify the computational results, the theoretical inversion of matrix (56) can be obtained as(57)Ytk=A−1tk=sintkcostk−costksintk∈ℝ2×2,which is given for checking the effectiveness and correctness of the three models and corresponding algorithms. The initial values of the example are arbitrarily set as(58)Y0=1.51.92.35.8,for all algorithms.


### 6.1.1. New Models of Zhang Neural Network

Firstly, we check the effectiveness and the correctness of ZMMMI model (10) with five discrete algorithms, i.e., ZMMMI2i algorithm (25), ZMMMI4i algorithm (27), ZMMMI6i algorithm (29), ZMMMI8i algorithm (31), and ZMMMI10i algorithm (33). The step size and task duration are uniformly set as h=0.005 and td=30 s. The corresponding results of computer experiments are shown in Figures 1 and 2. Figure 1(a) shows the residual errors Xk+1Ak+1−IF of the five ZMMMI algorithms, with ϵ=0.001 s and h=0.005. Figure 1(b) illustrates the residual errors Xk+1−Ak+1−1F of the same five ZMMMI algorithms as before. From the trajectories of all entries in Figures 2(a)–2(d), the solution of the model coincides with the theoretical solution perfectly. In addition, we see that the residual error trajectories synthesized by different ZMMMI algorithms quickly stabilized to the steady-state error level, after undergoing the initial hundreds of recursions.Figure 1
Residual error trajectories synthesized by five ZMMMI algorithms for example 1 (withϵ=0.001 s and γ=5.0). (a) Trajectories of Xk+1Ak+1−IF. (b) Trajectories of Xk+1−Ak+1−1F.
(a)(b)Figure 2
Solution trajectories synthesized byYk+1 (i.e., Ak+1−1) and by ZMMMI10i algorithm (33) for example 1 (with ϵ=0.001 s and γ=5.0). (a) Trajectories of y11 and x11. (b) Trajectories of y12 and x12. (c) Trajectories of y21 and x21. (d) Trajectories of y22 and x22.
(a)(b)(c)(d)


### 6.1.2. PGMDS Model

First, we check the effectiveness and the correctness of PGMDS model (38) with five discrete algorithms, i.e., PGMDS2i algorithm (40), PGMDS4i algorithm (41), PGMDS6i algorithm (42), PGMDS8i algorithm (43), and PGMDS10i algorithm (44). The step size and task duration are uniformly set as h=0.005 and td=30 s. The corresponding results of computer experiments are shown in Figures 3 and 4. Figure 3(a) shows the residual errors Xk+1Ak+1−IF of the five PGMDS algorithms, with ϵ=0.001 s and h=0.005. Figure 3(b) illustrates the residual errors Xk+1−Ak+1−1F of the same five PGMDS algorithms as before. From the trajectories of all entries in Figures 4(a)–4(d), the solution of the model coincides with the theoretical solution perfectly. In addition, we see that the residual error trajectories synthesized by different PGMDS algorithms quickly stabilized to the theoretical error level, after undergoing the initial hundreds of recursions.Figure 3
Residual error trajectories synthesized by five PGMDS algorithms for example 1 (withϵ=0.001 s and γ=5.0). (a) Trajectories of Xk+1Ak+1−IF. (b) Trajectories of Xk+1−Ak+1−1F.
(a)(b)Figure 4
Solution trajectories synthesized byYk+1 (i.e., Ak+1−1) and by PGMDS10i algorithm (39) for example 1 (with ϵ=0.001 s and γ=5.0). (a) Trajectories of y11 and x11. (b) Trajectories of y12 and x12. (c) Trajectories of y21 and x21. (d) Trajectories of y22 and x22.
(a)(b)(c)(d)


### 6.1.3. DDD Model

First, we check the effectiveness and correctness of DDD model (49) with five discrete algorithms, i.e., DDD2i algorithm (51), DDD4i algorithm (52), DDD6i algorithm (53), DDD8i algorithm (54), and DDD10i algorithm (55). The step size and task duration are uniformly set as h=0.0001 and td=30 s. The corresponding results of computer experiments are shown in Figures 5 and 6. Figure 5(a) shows the residual errors Xk+1Ak+1−IF of the five DDD algorithms, with ϵ=0.001 s and h=0.005. Figure 5(b) illustrates the residual errors Xk+1−Ak+1−1F of the same five DDD algorithms as before. From the trajectories of all entries in Figures 6(a) and 6(b), even for the model with the highest accuracy, DDD10i algorithm (55), its solution still cannot converge to the theoretical solution.Figure 5
Residual error trajectories synthesized by five DDD algorithms for example 1 (withϵ=0.001 s and γ=5.0). (a) Trajectories of Xk+1Ak+1−IF. (b) Trajectories of Xk+1−Ak+1−1F.
(a)(b)Figure 6
Solution trajectories synthesized byYk+1 (i.e., Ak+1−1) and by DDD10i algorithm (55) for example 1 (with ϵ=0.001 s and γ=5.0). (a) Trajectories of y11 and x11. (b) Trajectories of y12 and x12. (c) Trajectories of y21 and x21. (d) Trajectories of y22 and x22.
(a)(b)(c)(d)


## 6.1.1. New Models of Zhang Neural Network

Firstly, we check the effectiveness and the correctness of ZMMMI model (10) with five discrete algorithms, i.e., ZMMMI2i algorithm (25), ZMMMI4i algorithm (27), ZMMMI6i algorithm (29), ZMMMI8i algorithm (31), and ZMMMI10i algorithm (33). The step size and task duration are uniformly set as h=0.005 and td=30 s. The corresponding results of computer experiments are shown in Figures 1 and 2. Figure 1(a) shows the residual errors Xk+1Ak+1−IF of the five ZMMMI algorithms, with ϵ=0.001 s and h=0.005. Figure 1(b) illustrates the residual errors Xk+1−Ak+1−1F of the same five ZMMMI algorithms as before. From the trajectories of all entries in Figures 2(a)–2(d), the solution of the model coincides with the theoretical solution perfectly. In addition, we see that the residual error trajectories synthesized by different ZMMMI algorithms quickly stabilized to the steady-state error level, after undergoing the initial hundreds of recursions.Figure 1
Residual error trajectories synthesized by five ZMMMI algorithms for example 1 (withϵ=0.001 s and γ=5.0). (a) Trajectories of Xk+1Ak+1−IF. (b) Trajectories of Xk+1−Ak+1−1F.
(a)(b)Figure 2
Solution trajectories synthesized byYk+1 (i.e., Ak+1−1) and by ZMMMI10i algorithm (33) for example 1 (with ϵ=0.001 s and γ=5.0). (a) Trajectories of y11 and x11. (b) Trajectories of y12 and x12. (c) Trajectories of y21 and x21. (d) Trajectories of y22 and x22.
(a)(b)(c)(d)


## 6.1.2. PGMDS Model

First, we check the effectiveness and the correctness of PGMDS model (38) with five discrete algorithms, i.e., PGMDS2i algorithm (40), PGMDS4i algorithm (41), PGMDS6i algorithm (42), PGMDS8i algorithm (43), and PGMDS10i algorithm (44). The step size and task duration are uniformly set as h=0.005 and td=30 s. The corresponding results of computer experiments are shown in Figures 3 and 4. Figure 3(a) shows the residual errors Xk+1Ak+1−IF of the five PGMDS algorithms, with ϵ=0.001 s and h=0.005. Figure 3(b) illustrates the residual errors Xk+1−Ak+1−1F of the same five PGMDS algorithms as before. From the trajectories of all entries in Figures 4(a)–4(d), the solution of the model coincides with the theoretical solution perfectly. In addition, we see that the residual error trajectories synthesized by different PGMDS algorithms quickly stabilized to the theoretical error level, after undergoing the initial hundreds of recursions.Figure 3
Residual error trajectories synthesized by five PGMDS algorithms for example 1 (withϵ=0.001 s and γ=5.0). (a) Trajectories of Xk+1Ak+1−IF. (b) Trajectories of Xk+1−Ak+1−1F.
(a)(b)Figure 4
Solution trajectories synthesized byYk+1 (i.e., Ak+1−1) and by PGMDS10i algorithm (39) for example 1 (with ϵ=0.001 s and γ=5.0). (a) Trajectories of y11 and x11. (b) Trajectories of y12 and x12. (c) Trajectories of y21 and x21. (d) Trajectories of y22 and x22.
(a)(b)(c)(d)


## 6.1.3. DDD Model

First, we check the effectiveness and correctness of DDD model (49) with five discrete algorithms, i.e., DDD2i algorithm (51), DDD4i algorithm (52), DDD6i algorithm (53), DDD8i algorithm (54), and DDD10i algorithm (55). The step size and task duration are uniformly set as h=0.0001 and td=30 s. The corresponding results of computer experiments are shown in Figures 5 and 6. Figure 5(a) shows the residual errors Xk+1Ak+1−IF of the five DDD algorithms, with ϵ=0.001 s and h=0.005. Figure 5(b) illustrates the residual errors Xk+1−Ak+1−1F of the same five DDD algorithms as before. From the trajectories of all entries in Figures 6(a) and 6(b), even for the model with the highest accuracy, DDD10i algorithm (55), its solution still cannot converge to the theoretical solution.Figure 5
Residual error trajectories synthesized by five DDD algorithms for example 1 (withϵ=0.001 s and γ=5.0). (a) Trajectories of Xk+1Ak+1−IF. (b) Trajectories of Xk+1−Ak+1−1F.
(a)(b)Figure 6
Solution trajectories synthesized byYk+1 (i.e., Ak+1−1) and by DDD10i algorithm (55) for example 1 (with ϵ=0.001 s and γ=5.0). (a) Trajectories of y11 and x11. (b) Trajectories of y12 and x12. (c) Trajectories of y21 and x21. (d) Trajectories of y22 and x22.
(a)(b)(c)(d)


## 6.2. Example 2

The second time-variant matrix is a3×3 real matrix which is shown as follows:(59)Atk=5+costkcostk−sintkcostk3+sintksintkcostk−costk2+sintk,where Atk∈ℝ3×3. To verify the computational results, we utilize the theoretical inverse of Atk, which is denoted as Ytk=yijtk. Because Ytk is too complicated, it is omitted here. The initial values of the example are arbitrarily set as(60)Y0=9.50.70.6−1.1121.3−1.60.29.0.The task duration is uniformly set astd=30 s.


### 6.2.1. ZMMMI Model

The results of the numerical experiments for ZMMMI2i algorithm (25), ZMMMI4i algorithm (27), ZMMMI6i algorithm (29), ZMMMI8i algorithm (31), and ZMMMI10i algorithm (33) are shown in Figures 7 and 8. The experimental results are satisfactory as expected.Figure 7
Residual error trajectories synthesized by five ZMMMI algorithms for example 2 (withϵ=0.0001 s and γ=2.0). (a) Trajectories of Xk+1Ak+1−IF. (b) Trajectories of Xk+1−Ak+1−1F.
(a)(b)Figure 8
Solution trajectories synthesized byYk+1 (i.e., Ak+1−1) and by ZMMMI10i algorithm (33) for example 2 (with ϵ=0.0001 s and γ=2.0). (a) Trajectories of y11 and x11. (b) Trajectories of y12 and x12. (c) Trajectories of y13 and x13. (d) Trajectories of y21 and x21. (e) Trajectories of y22 and x22. (f) Trajectories of y23 and x23. (g) Trajectories of y31 and x31. (h) Trajectories of y32 and x32. (i) Trajectories of y33 and x33.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


### 6.2.2. PGMDS Model

The results of the numerical experiments for PGMDS2i algorithm (40), PGMDS4i algorithm (41), PGMDS6i algorithm (42), PGMDS8i algorithm (43), and PGMDS10i algorithm (44) are illustrated in Figures 9 and 10. The experimental results are satisfactory as expected.Figure 9
Residual error trajectories synthesized by five PGMDS algorithms for example 2 (withϵ=0.001 s and γ=5.0). (a) Trajectories of Xk+1Ak+1−IF. (b) Trajectories of Xk+1−Ak+1−1F.
(a)(b)Figure 10
Solution trajectories synthesized byYk+1 (i.e., Ak+1−1) and by PGMDS10i algorithm (44) for example 2 (with ϵ=0.001 s and γ=5.0). (a) Trajectories of y11 and x11. (b) Trajectories of y12 and x12. (c) Trajectories of y13 and x13. (d) Trajectories of y21 and x21. (e) Trajectories of y22 and x22. (f) Trajectories of y23 and x23. (g) Trajectories of y31 and x31. (h) Trajectories of y32 and x32. (i) Trajectories of y33 and x33.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


### 6.2.3. DDD Model

The results of the numerical experiments for DDD2i algorithm (51), DDD4i algorithm (52), DDD6i algorithm (53), DDD8i algorithm (54), and DDD10i algorithm (55) are shown in Figures 11 and 12. In terms of convergence and accuracy, the experimental results in this example are consistent with those in example 1.Figure 11
Residual error trajectories synthesized by five DDD algorithms for example 2 (withϵ=0.001 s and γ=1.0). (a) Trajectories of Xk+1Ak+1−IF. (b) Trajectories of Xk+1−Ak+1−1F.
(a)(b)Figure 12
Solution trajectories synthesized byYk+1 (i.e., Ak+1−1) and by DDD10i algorithm (55) for example 2 (with ϵ=0.001 s and γ=1.0). (a) Trajectories of y11 and x11. (b) Trajectories of y12 and x12. (c) Trajectories of y13 and x13. (d) Trajectories of y21 and x21. (e) Trajectories of y22 and x22. (f) Trajectories of y23 and x23. (g) Trajectories of y31 and x31. (h) Trajectories of y32 and x32. (i) Trajectories of y33 and x33.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


## 6.2.1. ZMMMI Model

The results of the numerical experiments for ZMMMI2i algorithm (25), ZMMMI4i algorithm (27), ZMMMI6i algorithm (29), ZMMMI8i algorithm (31), and ZMMMI10i algorithm (33) are shown in Figures 7 and 8. The experimental results are satisfactory as expected.Figure 7
Residual error trajectories synthesized by five ZMMMI algorithms for example 2 (withϵ=0.0001 s and γ=2.0). (a) Trajectories of Xk+1Ak+1−IF. (b) Trajectories of Xk+1−Ak+1−1F.
(a)(b)Figure 8
Solution trajectories synthesized byYk+1 (i.e., Ak+1−1) and by ZMMMI10i algorithm (33) for example 2 (with ϵ=0.0001 s and γ=2.0). (a) Trajectories of y11 and x11. (b) Trajectories of y12 and x12. (c) Trajectories of y13 and x13. (d) Trajectories of y21 and x21. (e) Trajectories of y22 and x22. (f) Trajectories of y23 and x23. (g) Trajectories of y31 and x31. (h) Trajectories of y32 and x32. (i) Trajectories of y33 and x33.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


## 6.2.2. PGMDS Model

The results of the numerical experiments for PGMDS2i algorithm (40), PGMDS4i algorithm (41), PGMDS6i algorithm (42), PGMDS8i algorithm (43), and PGMDS10i algorithm (44) are illustrated in Figures 9 and 10. The experimental results are satisfactory as expected.Figure 9
Residual error trajectories synthesized by five PGMDS algorithms for example 2 (withϵ=0.001 s and γ=5.0). (a) Trajectories of Xk+1Ak+1−IF. (b) Trajectories of Xk+1−Ak+1−1F.
(a)(b)Figure 10
Solution trajectories synthesized byYk+1 (i.e., Ak+1−1) and by PGMDS10i algorithm (44) for example 2 (with ϵ=0.001 s and γ=5.0). (a) Trajectories of y11 and x11. (b) Trajectories of y12 and x12. (c) Trajectories of y13 and x13. (d) Trajectories of y21 and x21. (e) Trajectories of y22 and x22. (f) Trajectories of y23 and x23. (g) Trajectories of y31 and x31. (h) Trajectories of y32 and x32. (i) Trajectories of y33 and x33.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


## 6.2.3. DDD Model

The results of the numerical experiments for DDD2i algorithm (51), DDD4i algorithm (52), DDD6i algorithm (53), DDD8i algorithm (54), and DDD10i algorithm (55) are shown in Figures 11 and 12. In terms of convergence and accuracy, the experimental results in this example are consistent with those in example 1.Figure 11
Residual error trajectories synthesized by five DDD algorithms for example 2 (withϵ=0.001 s and γ=1.0). (a) Trajectories of Xk+1Ak+1−IF. (b) Trajectories of Xk+1−Ak+1−1F.
(a)(b)Figure 12
Solution trajectories synthesized byYk+1 (i.e., Ak+1−1) and by DDD10i algorithm (55) for example 2 (with ϵ=0.001 s and γ=1.0). (a) Trajectories of y11 and x11. (b) Trajectories of y12 and x12. (c) Trajectories of y13 and x13. (d) Trajectories of y21 and x21. (e) Trajectories of y22 and x22. (f) Trajectories of y23 and x23. (g) Trajectories of y31 and x31. (h) Trajectories of y32 and x32. (i) Trajectories of y33 and x33.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


## 6.3. Example 3

This example compares solely the effective models, i.e., ZMMMI (10) and PGMDS (38), with a 4×4 real matrix, which is shown as follows:(61)Atk=6+ckckckckck5+sk−skckckck4+sk−skskskck2+sk,where Atk∈ℝ4×4 and sk and ck denote sintk and costk, respectively. To verify the computational results, we utilize the theoretical inverse of Atk, which is denoted as Ytk=yijtk. Because Ytk is too complicated, it is omitted here. The initial values of the example are arbitrarily set as(62)Y0=6111150111400012.The task duration is uniformly set astd=30 s.


### 6.3.1. ZMMMI Model

The results of the numerical experiments for ZMMMI2i algorithm (25), ZMMMI4i algorithm (27), ZMMMI6i algorithm (29), ZMMMI8i algorithm (31), and ZMMMI10i algorithm (33) are shown in Figures 13(a) and 13(b). The experimental results are satisfactory as expected.Figure 13
Residual error trajectories synthesized by five ZMMMI algorithms for example 3 (withϵ=0.0001 s and γ=2.0). (a) Trajectories of Xk+1Ak+1−IF. (b) Trajectories of Xk+1−Ak+1−1F.
(a)(b)


### 6.3.2. PGMDS Model

The results of the numerical experiments for PGMDS2i algorithm (40), PGMDS4i algorithm (41), PGMDS6i algorithm (42), PGMDS8i algorithm (43), and PGMDS10i algorithm (44) are illustrated in Figures 14(a) and 14(b). The experimental results are satisfactory as expected.Figure 14
Residual error trajectories synthesized by five PGMDS algorithms for example 3 (withϵ=0.002 s and γ=5.0). (a) Trajectories of Xk+1Ak+1−IF. (b) Trajectories of Xk+1−Ak+1−1F.
(a)(b)


### 6.3.3. Comparison and Discussion

We mainly investigate the residual error of the models and the coincidence between the solution matrix with the ground-truth matrix inverse. From the discrete simulation results of two examples of each five discrete algorithms, it can be seen that the convergence of ZMMMI model (10) and PGMDS model (38) is good, which is completely consistent with the conclusion of Proposition 1 and Theorems 1 and 2. Therefore, it can be concluded that the two models and corresponding discrete algorithms are effective. However, in the case of DDD model (49), the residual errors and the coincidence of the solution matrix entries are not satisfactory, so we evidently summarize that the effectiveness of DDD model (49) and corresponding discrete algorithms is not enough.


### 6.3.4. Remark

Note that the convergence of DDD models is shown in [12, 47], where the DDD models are utilized to solve the time-varying nonlinear optimization problems. There are three points worth further discussing. First, in [12], the initial value is set to be the theoretical value, and the computational process is not disturbed. Second, comparatively, in this study, the initial value is not close enough to the theoretical value. Third, the higher-order discrete algorithm (i.e., the ten-instant discrete algorithm) is more easily affected by the rounding error disturbance [44]. Thus, the experimental results of the DDD model in this study showing the divergence are actually complementing the previous research studies, in addition to the confirmation of [47] about divergence. Therefore, we summarize that the DDD model is generally less effective, with further in-depth investigation being also a future research direction.


## 6.3.1. ZMMMI Model

The results of the numerical experiments for ZMMMI2i algorithm (25), ZMMMI4i algorithm (27), ZMMMI6i algorithm (29), ZMMMI8i algorithm (31), and ZMMMI10i algorithm (33) are shown in Figures 13(a) and 13(b). The experimental results are satisfactory as expected.Figure 13
Residual error trajectories synthesized by five ZMMMI algorithms for example 3 (withϵ=0.0001 s and γ=2.0). (a) Trajectories of Xk+1Ak+1−IF. (b) Trajectories of Xk+1−Ak+1−1F.
(a)(b)


## 6.3.2. PGMDS Model

The results of the numerical experiments for PGMDS2i algorithm (40), PGMDS4i algorithm (41), PGMDS6i algorithm (42), PGMDS8i algorithm (43), and PGMDS10i algorithm (44) are illustrated in Figures 14(a) and 14(b). The experimental results are satisfactory as expected.Figure 14
Residual error trajectories synthesized by five PGMDS algorithms for example 3 (withϵ=0.002 s and γ=5.0). (a) Trajectories of Xk+1Ak+1−IF. (b) Trajectories of Xk+1−Ak+1−1F.
(a)(b)


## 6.3.3. Comparison and Discussion

We mainly investigate the residual error of the models and the coincidence between the solution matrix with the ground-truth matrix inverse. From the discrete simulation results of two examples of each five discrete algorithms, it can be seen that the convergence of ZMMMI model (10) and PGMDS model (38) is good, which is completely consistent with the conclusion of Proposition 1 and Theorems 1 and 2. Therefore, it can be concluded that the two models and corresponding discrete algorithms are effective. However, in the case of DDD model (49), the residual errors and the coincidence of the solution matrix entries are not satisfactory, so we evidently summarize that the effectiveness of DDD model (49) and corresponding discrete algorithms is not enough.


## 6.3.4. Remark

Note that the convergence of DDD models is shown in [12, 47], where the DDD models are utilized to solve the time-varying nonlinear optimization problems. There are three points worth further discussing. First, in [12], the initial value is set to be the theoretical value, and the computational process is not disturbed. Second, comparatively, in this study, the initial value is not close enough to the theoretical value. Third, the higher-order discrete algorithm (i.e., the ten-instant discrete algorithm) is more easily affected by the rounding error disturbance [44]. Thus, the experimental results of the DDD model in this study showing the divergence are actually complementing the previous research studies, in addition to the confirmation of [47] about divergence. Therefore, we summarize that the DDD model is generally less effective, with further in-depth investigation being also a future research direction.


## 7. Conclusion

In this paper, we have shed some light on the matrix inversion solution models derivation, i.e., ZMMMI model (10), PGMDS model (38), and DDD model (49), from PEs. It has provided a new perspective to make full use of the theoretical value of PEs. First, with the substitution technique and design formula of ZNN, we have investigated and proposed the new model of ZNN for matrix inversion problem. Then, we have discussed the convergence and accuracy of ZMMMI model (10) and presented two theorems and proofs about it. On the basis of the model, we have developed five ZMMMI algorithms to discretize continuous-time ZMMMI model (10). Second, with the substitution technique and design formula of ZNN, we have investigated and presented PGMDS model (38) for matrix inversion problem. On the basis of the model, we have shown five PGMDS algorithms to discretize continuous-time PGMDS model (38). Third, with the substitution technique and design formula of ZNN, we have presented DDD model (49) for matrix inversion problem. On the basis of the model, we have developed five DDD algorithms to discretize continuous-time DDD model (38). Fourth, we have prepared three examples to calculate the inverse of matrices by using the above three models, respectively. The results illustrate that ZMMMI model (10) and PGMDS model (38) are effective, while DDD model (49) is less effective. In the future, based on PEs, we will use the substitution technique and ZNN design formula to discuss more efficient RNN models for TVMI problem, including pseudo-inversion and complex matrix inversion.



---

*Source: 4227512-2021-11-20.xml*