**5. Illustrative Examples**

In this section we consider some numerical examples to demonstrate the performance of our RNN implementation for solving a system of linear algebraic equations.

#### *5.1. Illustrative Example 1*

Let us define M and V as follows: Here, instead of using the identity matrix in Equation (1) we use a vector V. This does result in a linear system of equations.

$$M = \begin{bmatrix} -3.0755 & 0.5004 & 3.8135 \\ 0.0739 & -2.2382 & -4.7532 \\ -4.7576 & 1.9624 & -1.5879 \end{bmatrix}, V = \begin{bmatrix} 7.920 \\ 7.983 \\ 9.633 \end{bmatrix} \tag{37}$$

Thus, one can find *X* as following:

$$X = M^{-1}V = \begin{bmatrix} -3.330 \\ -3.305 \\ -0.175 \end{bmatrix} \tag{38}$$

This above given value of X is the target value (ground truth). In the following, we shall see how far how fast the models developed (Equation (24) and/or Equation (31)) do reach well this target value.

M is not singular and therefore can be inverted. Thus, our system of equations has one unique solution. C, as explained in the previous section, is a weight values-matrix which is calculated using the following formula:

$$C = \sum\_{i=0}^{n} \left(\boldsymbol{M}^T \boldsymbol{M}\right)^i \tag{39}$$

Increasing the value of n will increase the convergence rate of Equation (24) or Equation (31) or (29). Here, we choose the value n = 3. The corresponding RNN parameters are defined as follows (see Figure 3):

$$\begin{aligned} \text{C} &= M^T M M^T M M^T M + M^T M M^T M + M^T M + I\\ \text{C} &= 10,000 \times \begin{bmatrix} 4.6283 & -2.2035 & -2.7439\\ -2.2035 & 1.3283 & 2.5763\\ -2.7439 & 2.5763 & 7.5178 \end{bmatrix} \\ &\quad X\_0 = \begin{bmatrix} 0\\ 0\\ 0 \end{bmatrix} \end{aligned} \tag{40}$$

These parameters are implemented in the previously described model in Simulink (see Figure 3) for numerical simulations and F (for simplicity) is taken as a linear function (Figure 1).

Figure 4 is showing the system simulation results during the time interval *t* = [0, 0.5]. The output (value of state vector) at the end of this time interval is the following:

$$X(0.5) = \begin{bmatrix} -3.329 \\ -3.304 \\ -0.1749 \end{bmatrix} \tag{41}$$

This above obtained result shows that the output of our system model is very close to the analytical solution and the difference is small and approximately 0.001, which can furthermore be corrected but adjusting the step size.

**Figure 4.** Result of our model simulation for solving the illustrative example with respect to solving an algebraic system of linear equations. Please note in this case we have taken n = 3 in Equation (24) or Equation (31).

#### *5.2. Illustrative Example 2*

Now we do the same experience for calculating the inverse of a time-varying matrix *M*. The matrix *M* is defined as follows:

$$M = \begin{bmatrix} \operatorname{Sin}(t) & -\operatorname{Cos}(t) \\ \operatorname{Cos}(t) & \operatorname{Sin}(t) \end{bmatrix} \tag{42}$$

the matrix *M* is always invertible and it does satisfy the conditions of Equation (1). Therefore, we can use it for testing our dynamic system model.

$$M^T = M^{-1} = \begin{bmatrix} \operatorname{Sin}(t) & \operatorname{Cos}(t) \\ -\operatorname{Cos}(t) & \operatorname{Sin}(t) \end{bmatrix} \tag{43}$$

Also, we define F as a linear function (Figure 1).

As we have explained previously (see Table 2), we calculate the weight *C* again for n = 3 as the following:

$$\mathbf{C} = \mathbf{M}^{\mathrm{T}} \mathbf{M} \mathbf{M}^{\mathrm{T}} \mathbf{M} \mathbf{M}^{\mathrm{T}} \mathbf{M} + \mathbf{M}^{\mathrm{T}} \mathbf{M} \mathbf{M}^{\mathrm{T}} \mathbf{M} + \mathbf{M}^{\mathrm{T}} \mathbf{M} + \mathbf{I} \tag{44}$$

Now all parameters are introduced in the Simulink model to simulate our system model (see Equation (24)). The simulation is done for *t* = [0, 2.0]. Figure 5 does contain benchmarking results for the first element of the time-varying matrix (i.e., the matrix element in first row, first column). Here (see Figure 5) the speed of convergence is compared to that of the original Zhang model. Thereby, for our model we consider di fferent values of the parameter n. One can clearly see that our novel model strongly outperforms the original Zhang model. Table 3 shows the di fference very clearly: Here, it is clear that by choosing n = 3 our model will converge to the problem solution 4 times faster than the ZNN model; notice that the ZNN model corresponds to n = 0 in Table 3. Amount of memory for saving this model is very small as we need to save main templates parameter of dynamic model in Equation (24) or Equation (31). Therefore, as see in Table 3, the memory usage is very small and it can be implemented on small computers.

**Table 3.** Comparison of performance of time varying matrix inversion with change of parameter n of Equation (24) or Equation (31). Please note n = 0 is ZNN model.


#### *5.3. Illustrative Example 3*

In the last illustrative experiment, we try to test and analyze the e ffect of changing the function F in Equation (31) on the convergence rate of our model.

For this test, we take the same parameter setting as in the first experience (see illustrative example 1), but we change the function F and calculate coe fficients for *n* = 3. The selected functions for this test are the following ones: linear function ( *Y* = *X*), sigmoid function ( *Y* = 1 2 (|*X* + 1| − |*X* − 1|), arctangent function ( *Y* = arctan(*X*)), tangent hyperbolic function ( *Y* = tan h(*X*)), and polynomial function ( *Y* = <sup>X</sup><sup>3</sup>). All other parameters are fixed such to show the correct properties of the functions. As we can see in the results of this simulation experiment (see Figure 6) the polynomial function is showing the very best convergence amongs<sup>t</sup> all considered functions. Therefore, it is evident that selecting the appropriate function can help our model to converge faster to the solution of the matrix inversion problem. It is also evident that higher values of n will perform much better.

Table 4 is showing our simulation result until t = 0.05 for this illustrative example 3. It does clearly show the e ffect of using di fferent functions on the system/model convergence. Both polynomial and linear functions are the best functions for this task, while both sigmoid and tanh functions are the worst functions to be used in this context for matrix inversion.

**Table 4.** Performance comparison of different function types used in Equation (31), whereby all are monotonic increasing functions. The polynomial function shows the better convergence compared to the other functions.

**Figure 5.** Showing the difference in convergence speed between the original Zhang method and our new developed method (see Equation (24)). The solid lines are the one obtained from the model solving and the dashed lines are the target values. (**a**) convergence of the original Zhang model; (**b**) until (**d**): convergence profiles of our novel method for different values of n. All graphs (a) until (d) are plotted for the first element of the time-varying matrix M (i.e., first row, first column).

**Figure 6.** Showing differences in convergence speed while considering different function types for F in Equation (31). Please note that n is 3 for the polynomial function F.

The superiority with respect to convergence speed of the polynomial function F is very evident.

#### **6. Comparison of Our Novel Method with Previous Studies**

In this section we compare our novel model with the following related methods which are well-known from the relevant literature: gradient descent method, Zhang dynamics, and Chen dynamics. Hereby, 1000 different static random matrices are generated for use in the experiment. Finally, we then sum up the results obtained as shown in Figure 7. Figure 7 does display how the error converges towards zero when the different models are executed over time; the intention is to hereby illustrate the convergence speed of all considered models. Hereby, the error function is calculated by using the formula *MX* − *I*.

By comparing the different methods with our method (Figure 7), the fastest convergence of our method to the exact solution is observed. Specifically, our method is at least 3 times faster than the Chen method while implemented on CPU (the results of Figure 7 were obtained on CPU). An implementation on a multiple-core platform should display a much higher relative and comparative speed-up performance of our model. Obviously, by using more coefficients (i.e., higher values of n) in Equation (24) we can reach a much better convergence rate.

**Figure 7.** Comparison of different DNN models with respect to convergence speed under the same conditions—benchmarking. In this figure, GD refers to the "gradient descent" method. Our novel method (see novel in the figure; Equation (24)) is calculated for n = 3. The fast convergence of our method to the exact solution is clearly demonstrated. It is also clear that this convergence would be much faster if higher values of n are taken (n > 3).
