Applications of the Network Simulation Method to Differential Equations with Singularities and Chaotic Behaviour

: In this paper, we deal with some applications of the network simulation method (NMS) to the non-linear differential equations derived of a parametric family associated to stated problems by Newton in and others like the parabolic mirror and van der Pol non-linear equation. We underly the efﬁcientcy of the (NMS) method, compare it with Matlab procedures and present ﬁgures of solutions of the equations obtained by it on the mentioned problems. Additionally, we introduce also the electric-electronic circuits we have designed to be able of obtaining the solutions of the referred equations.


Introduction
In many problems in science and engineering the network simulation method NSM [1] has demonstrated to be very efficient, for instance, in the simulation of the Frenkel-Kontorova-Tolinson model on friction of two sliding bodies [2].
Friction is a phenomenon present in all mechanical devices in which relative displacements happen. There are different cases to be considered. One of them is the dry friction. In this case the main magnitude is the force necessary to keep a sliding body moving at a given velocity from another considered motionless or velocity threshold at which this force is null.
Solutions in the dynamical problem are very sensitive on the parameters involved in the problems. Such parameters are the interaction force amplitude between atoms of the bodies, the stiffness between the atoms and bulk velocity.
Such solutions are a good example of phenomenon where they can be obtained of an approximate way using NSM with a numerical process shorter than other traditional methods offered in the literature.
In this paper we apply the NSM for obtaining numerical solutions of three examples, a parametric family of non-linear differential equations, the problem of the parabolic mirror and the well-known van der Pol equation. Such approach addresses several difficulties like singularities or chaotic behaviours. The solutions are compared with the results obtained by Matlab. Additionally, the NSM can be implemented in a very efficient and open software used in electronic circuit design, called Ngspice, which has the advantage of requiring a short computing times in comparison with other mathematical software.
The design of the network model implies the following steps: • Select the equivalence between physical and electrical variables; • Apply the initial and boundary conditions using the appropriate electrical devices. The open-source free-software called Ngspice, used in circuit design as numerical simulation software, allows us short computing time by a suitable time step from the point of view of convergence which is adjusted continuously.

Methodology
The selected problems are a family of differential equations with infinite slope in some points, the differential equation of a parabolic mirror and the differential equation for the behaviour of van der Pol oscillator.
These equations whose solutions are examples of functions setting out singularities and chaotic behaviour and will be useful to check the capacity of the NSM.
The network simulation method consists in transforming the studied system into electrical circuits with the same governing equations. Physical variables from original system have their equivalent in electrical variables. For example, the variable y in the family of differential equations with infinite slope in some points is equivalent to the voltage in a condenser node, as it will be explained more clearly in the applications of the NSM. Once the circuits equations are solved, using the traditional numerical integration method but with the same efficient techniques used with very-large-scale integration (VLSI), we will have the solution of the physical problem. Only it is necessary to consider again the equivalence between the electrical and the physical variables. Thus, the NGSpice software allows us to simulate electric circuits and analyse the response of them.
In the works of Solano et al. [3] there is a more detailed description of the precision of the NSM in some cases where it is difficult to get the convergence. The equivalent electrical network used for this method is solved by Kirchoff's laws and a convenient time step considering a stable convergence [4][5][6][7]. In the aforementioned paper, the interaction between the parameters called relative tolerance and chaotic behaviour of the system is analysed. As results of that study the solution stability is guaranteed.
The algorithm used by the software reaches the convergence for any iteration k when the following condition is achieved: where v n is the voltage in the node n and V n max = MAX |v k+1 n |, |v k n | These two parameter, RELTOL and VNTOL, are related to the time step. The values of RELTOL and VNTOL is chosen by the user considering his experience from other solved problems or checking the same problem with several values of these parameters. Besides, several works have been devoted to truncation error using PSpice, from which NGSpice comes from [8][9][10].
From the Vladimeresku analysis using circuits with condenser in which the exact solution is available, Forward-Euler (FE) expression is more liable to the value of time step than the Backward-Euler (BE) expression. The criteria used by software as Ngspice is that the Local Truncation Error (LTE) of constitutive equations for each branch must be limited.
The FE expression provides an inappropriate solution when the time step exceeds a certain value related to the capacitance of the condenser, τ, used in the analysed circuit. The voltage in the nodes belonging to the analysed circuit in which the right solution is zero is assessed by BE expression as a value which decreases over time. On the other hand, the Trapezoidal (TR) procedure provided a solution which oscillates with a decreasing amplitude around zero if time step is greater than 2τ. The time constants associated to the electrical networks have time constants utterly different, thus their governing equations define stiff systems. The methods aforementioned must be stiffly stable, that means that the time step can be higher than any time constant of the electrical network. Only the explicit methods are not stiffly stable. The stability properties of the Gear method with order higher than 1 is good. The simulation provides a solutions according to the selected integration algorithm. The characteristics of Gear and TR algorithms are different despite of the coincidence in the solutions. Speaking about the TR algorithm convergence it was mentioned his oscillatory behaviour. This oscillation is related with the fact that the time step are not smaller than a threshold value. Conversely, the Gear expression of order 2 introduces damping in the solution.
Another advantage of NSM is that it supplies relationships between y (t) and y (t) in cases of second order differential equations. In most cases NSM gives the solutions is easier than other metohds like Matlab. In order to obtain the solution by numerical methods of a second-order differential equation like Runge-Kutta, the second-order equation must be transformed into a system of first-order equations with the difficulty that this implies.

A Family of Differential Equations
We are interested in the analysis of solutions of the components of the following family of differential equations depending on a non negative real parameter a. Some of them are formulations associated to problems stated by Newton on the movement of bodies in the space [11].
In the family, there is a common map F(t, y) = t/(1 − y) plus a parameter a. Such parameter increases the slopes of the curves which represent the solutions.
is a real analytic function defined in the open rectangle: and for |y(t)| < 1 it is t 1−y(t) = 1 + t(1 + y(t) + y(t) 2 + · · · + · · · + y(t) n + · · · ) = 1 + t ∑ ∞ n=0 y(t) n . When a = 0 the equation is the unique of the family which can be integrated by the separation variables method obtaining the unicity of solution for real analytic functions by application the Cauchy theorem on existence and unicity of solutions [12] y Figure 1 shows this solution where two methods NSM and Matlab, give the same graphic y(t) for t < 1. The circuits designed for NSM are represented in Figures 2 and 3.
When a > 0 it is not possible the explicit integration of such equations. However, the use of power series allows us to solve all cases. Such approach was introduced by Newton [13].  Proof. Since F a (t) = a + t 1−y(t) is real analytic in the rectangle set S a (t, y) = {(t, y) : 0 < t < t a , y < 1} using the theorem of Cauchy [9], for each a > 0 there is unique solution which is real analytic. Then: y a (t) = i=∞ i=0 a i t i where a i > 0 are real numbers and the power series is absolutely convergent. The coefficients a i can be obtained recursively if we have obtained y a (t), y a (t) 2 , ..., y a (t) n , ...   Proof. Since F a (t) = a + t 1−y(t) is real analytic in the rectangle set S a (t, y) = {(t, y) : 0 < t < t a , y < 1} using the theorem of Cauchy [9], for each a > 0 there is unique solution which is real analytic. Then: y a (t) = i=∞ i=0 a i t i where a i > 0 are real numbers and the power series is absolutely convergent. The coefficients a i can be obtained recursively if we have obtained y a (t), y a (t) 2 , ..., y a (t) n , ... Solutions for a > 0 are given in the following result where with y a (t) we denote the solution for each a: Proof. Since F a (t) = a + t 1−y(t) is real analytic in the rectangle set S a (t, y) = {(t, y) : 0 < t < t a , y < 1} using the theorem of Cauchy [12], for each a > 0 there is unique solution which is real analytic. Then: y a (t) = ∑ i=∞ i=0 a i t i where a i > 0 are real numbers and the power series is absolutely convergent.
The coefficients a i can be obtained recursively if we have obtained y a (t), y a (t) 2 , . . . , y a (t) n , . . . For this we are using the Cauchy product of series.
In fact, we have Mathematics 2021, 9, 1442 5 of 12 where p i = a 0 a i + a 1 a i−1 + . . . + a i a 0 and i ∈ N (8) and The recurrence to obtain the above coefficients come from the equality y a (t) = a 1 + 2a 2 t + 3a 3 t 2 + . . . + na n t n−1 + . . .
We have generated the coefficients automatically by using a recurrent method, Figure 4. It is not difficult to see that the coefficients are the components of a sequence of positive real numbers convergent to zero. Therefore, in the interval (0, y a (t a )), where t a is the moment in which y reaches the value 1, the power series y(t) exists and the same for all its derivatives. The variable y(t) is increasing in the former interval and since y (t) is positive the map is convex in the same interval. In t = 0 is y (0) = a and lim t→α − a (y(t)) = +∞. Figures 5 and 6 show the results from, Matlab, and the NSM by Ngspice. The solutions become different for t > 0.5.
For t > t a it is not possible to obtain an explicit solution of any equation. Therefore, we have to use numerical approximations to solve it and to study and interpret graphically their behaviour.
The circuits designed for NSM for this case are represented in Figures 7 and 8. A particular case of the family (3) is when a = 1 ; y(0) = 0 When t ≥ t a from Figure 5 we observe the oscillation around y = 1.        In the interval (t a , t 1 ), where t 1 , t 2 , . . . , t i are the times greater than t a in which the slope of the function y(t) tends to infinity, the curve is decreasing holding y a (t a ) = −∞ and y a (t 1 ) = +∞ and increasing in (t 1 , t 2 ) with y a (t 1 ) = +∞ and y a (t 2 ) = −∞. The shape of the solution in the interval (0.5, 0.9) for Matlab code can be seen in the Figure 9, and in the interval (0.5446, 0.5449) for NSM in the Figure 10. Additionally, in t a , t 1 , t 2 , . . . the curves are non-differentiable. Additionally, in t a , t 1 , t 2 , . . . the curves are non-differentiable.
After zooming in on both Figures 9 and 10 the branches obtained by the NSM showing the convexity changes more precise than using Matlab code.  Now we will deal with other properties of the solutions of the equations of the family of differential equations when t ≥ t a . In fact, we use the notion of total variation denoted by V J of a continuous map f : J ⊂ R where J is compact. In fact, we follow the definition in the work of Hewitt and Stromberg [14].
is a finite partition of J being x 0 and x n P the extremes of J.
There is a connection between the total variation of a map f and the property of sensitive dependence on initial conditions. It means that if we take > 0 and d = |y 0 − y 1 | < where y 0 , y 1 are initial conditions for t = 0, there exists a t such that for all t ≥ t the distance d is greater than . The following result can be seen in [15] for discrete dynamical systems on a compact interval of the real line. However, doing simple changes in the statement of the result we adapt the result to differential equations. In fact, we obtain Proposition 1. The family (3) depending on a has sensitive dependence on initial conditions and it holds lim t→∞ (y(t)) = ∞ (20) for every closed subinterval of [α, β] for α > t a .

Solution of Differential Equation of a Parabolic Mirror
In the well-known book by Newton [11], we can find the definition of parabolic mirror. By symmetry, such a Mirror will have in space R 3 the shape of a surface of revolution generated by rotating an APB curve as indicated in the Figure 11. We can applied the theorem on existence and uniqueness of solutions used in Theorem 1.

Solution of differential equation of a parabolic mirror 103
In the well-known book by Newton [1] we can find the definition of parabolic mirror.

104
By symmetry, such a Mirror will have in space R 3 the shape of a surface of revolution 105 generated by rotating an APB curve as indicated in the Figure 10.

106
The ray that exits with angle θ exits parallel to the axis OX as indicated in the Figure   107 10.
What is the homogeneous differential equation that verifies the APB curve. Rewriting the (22) in the following way, is ± x 2 + y 2 = x + c with c ∈ R and simplifying x 2 + y 2 = x 2 + 2xc + c 2 , it results that it is the equation of all the focus parabolas at the origin and axis OX. Now, we try to solve the (22) with the NSM. In this method, the spatial variable x is replaced by t, the time: 2ty + yy 2 − y = 0 (26) The software used only provides the solution for t ≥ 0. In order to obtain the solution for t ≤ 0 it is necessary to turn the curve around the y axis, which is equivalent to multiply t by negative sign. In this way, the (26) is written as What is the homogeneous differential equation that verifies the APB curve. Rewriting the (22) in the following way, is ± x 2 + y 2 = x + c with c ∈ R and simplifying x 2 + y 2 = x 2 + 2xc + c 2 , it results that it is the equation of all the focus parabolas at the origin and axis OX. Now, we try to solve the (22) with the NSM. In this method, the spatial variable x is replaced by t, the time: The software used only provides the solution for t ≥ 0. In order to obtain the solution for t ≤ 0 it is necessary to turn the curve around the y axis, which is equivalent to multiply t by negative sign. In this way, the (26) is written as If we turn around the Y axis, the solution from the corresponding circuit and call x to t, the Figure 14 is obtained.

Van der Pol Equation
The well-known differential equation of van der Pol oscillator without considering excitation force is: where a > 0. The equation for the same system considering sinusoidal excitation force is: where a = 8.53 is a typical value in the literature. For a detailed treatment of Equations (29) and (30) see [16]. We apply a theorem of existence and uniqueness to see that the equation fulfils the theorem. Now, only numerical approach to solution is possible to be done. That is why we directly apply the NSM method.
In this case, the designed network for (30) is shown in Figure 15 and the solution in Figures 16 and 17. For a > 1 the equation is stiff. In Figure 16 The well-known differential equation of Van der Pol oscillator without considering excitation force is: where a > 0. The equation for the same system considering sinusoidal excitation force is: In this case the designed network for (30) is shown in Figure 14 and the solution in

Conclusions
The former examples show the capacity of the NSM to address equations with some singularities or chaotic behaviour like the mentioned family of differential equations. Besides, the use of this method in case of system of equations can be implemented with more easiness as other methods because it is only necessary to copy the networks already designed with the only care to change the denomination of each network.
Additionally, it can be implemented also for solving non-linear systems of differential equations.