Convergence and stability criteria for numerical solutions of partial differential equations in science and engineering

. - This paper explores the computational aspects of simulation and modeling applied to the solution of the Heaviside equation, considering the relevance of the stability and convergence of the solutions. For this purpose, a second-order finite difference scheme was implemented as the primary approach for studying atmospheric discharges (ATDI). The programming language used was Matlab, which facilitated calculating the induced currents in the study scenario. The centered, forward, and backward finite difference approaches were considered for the numerical implementation. System validation tests were performed to demonstrate the effectiveness of the design and convergence to the second order with the centered difference approach.


I. INTRODUCTION
The stability and convergence of differential equation solutions represent two essential pillars in any simulation process.Unfortunately, on many occasions, the numerical solution stability does not receive the necessary attention or is only superficially addressed in specific experiments.However, this study focuses rigorously on these aspects of vital importance in the context of simulations.The analysis focuses on the propagation of waves related to the phenomenon of atmospheric discharges, a problem that directly impacts electric utilities.The influence of this phenomenology on the extensive network of transmission lines and devices in electrical substations causes interruptions in power distribution and results in economic losses due to sharp fluctuations in current and voltage magnitudes, significantly when they exceed the isolation thresholds of such devices.
The intricate nature of this phenomenon complicates obtaining analytical solutions for the equations that describe it, giving rise to the imperative need to resort to numerical models supported by increasingly advanced software and hardware.Various numerical methods have paved the way for solving several applications linked to atmospheric discharges, thus obtaining highly accurate numerical solutions.Extending previous approaches, this study delves into the computational aspects of significant relevance in the simulation process.In particular, special emphasis is given to convergence and stability tests, whose conclusive results ensure the reliability of the solutions obtained.This work sheds light on the critical importance of stability and convergence in simulations and contributes to a more adequate understanding and management of the issues associated with lightning strikes.By delving deeper into these essential aspects, a solid foundation is laid for future research and advances in mitigating the negative impacts of this natural electrical phenomenon.

II. THE PHYSICAL MODEL AND THE MATHEMATICAL MODEL
The system under analysis is a relationship between the atmosphere and the earth, represented as a cylindrical condenser.The separation between these entities is defined by a discharge channel, which is assumed straight and without deviations.The cylindrical condenser has a distance of 1 km between its components, while the discharge channel has a diameter of approximately 20 cm.During the simulation of this scenario, a characteristic altitude of 1000 m is considered.This altitude is divided into 250 segments, each with a time interval Δt of 0.5 µs, for a detailed representation of the process.The system of equations used in this work is the one proposed in references [1] and [2], from which (1) and ( 2) are extracted as main equations where V(x,t) is the cloud voltage, I(x,t) is the discharge current and R, L, C and G are the resistance of the channel, the inductance associated with the discharge, and the inductance associated with the discharge, the capacitance and G the conductance, respectively.Combining equations (1) and (2) yields the well-known telegrapher's equation for current: It is easy to show that equation (3) has the same form for the lightning voltage.In this case, it is solved numerically to determine the return discharge current with the appropriate initial and boundary conditions.In addition to the actual parameters that allow the atmospheric discharge to be accurately described, it should be noted that this equation has been solved on other occasions, and the authors make assumptions that greatly simplify it [6].

III. NUMERICAL DEVELOPMENT AND IMPLEMENTATION
For the discretization of equation ( 3), the basic idea of finite differences is to replace continuous space-time with a discrete set of points.The time step between two consecutive levels is denoted ∆t, and the distance between two adjacent points in space ∆x [3].The next step is to replace the differential equations with discrete equations.This is achieved by approximating the differential operators by their corresponding finite differences and considering the function values at adjacent points in the mesh.In this way, an algebraic equation is obtained at each point of the grid for each term of the differential equation.These equations involve the function values at the moment considered and their nearest neighbors.Using the approximations obtained by Taylor's series development, its corresponding discretization replaces each differential operator in equation ( 3), and the following is obtained.
This equation has the property that it involves only two values of the wave function at the last time level; the value , then, can be cleared in terms of values at previous times to obtain the equation of the discretized telegrapher, which will later be implemented in Matlab.
To numerically implement this equation, it is necessary to specify the boundary and initial condition.The current in the cloud is assumed to be a constant value taken as zero.Therefore, the cofounder's condition is where h is the height of the cloud where the discharge starts, at t = 0, the current distribution is assumed to be a constant value.Therefore, the initial condition is

B. Stability and convergence analysis
For any simulation, stability and convergence tests are necessary, among other aspects, because of the following [3,4,5]: When applying finite difference methods involving numerical solutions of second-order differential equations, in addition to considering errors due to rounding and series truncation, it is essential to remember that finite difference expressions are derived from a Taylor series approximation up to the second term.This approximation introduces an error in the calculations that need to be controlled.
There are also errors associated with the limitations of fitting the physical model to accurately describe the natural system and the constraints imposed by the boundary conditions.Convergence and stability validation verify the correctness of the approximations made by finite differences in the discretization process.
Methods are available to evaluate the convergence of a code when analytical solutions are known.In this study, the electrical parameters are considered to achieve a sufficiently refined mesh, thus allowing code convergence in as few iterations as possible.To achieve this, the CFL condition is applied to ensure that the mesh will enable us to obtain a solution closer to reality.
A particular case arises when time-stepping schemes are used as a numerical solution method.Consequently, in many simulations, the time step must be smaller than a characteristic value to prevent the code from becoming unstable and producing incorrect results.This constraint is commonly called the numerical propagation speed in the CFL condition.
To assess stability, it is sufficient to ensure that the solutions of the differential equations do not exhibit oscillations, which can be determined by observing the graphs representing the time evolutions.This requires allowing the time-dependent variables to evolve over a prolonged period to verify that they do not exhibit unexpected oscillations or fluctuations.
Convergence refers to how the numerical solution approaches a known or actual analytical solution.In situations where the analytical solution is unknown, it is possible to use an iterative approach to evaluate convergence.The numerical value of the first iteration is used as the analytical approximation and compared to the numerical value of the second iteration.Then, the value of the second iteration is taken as the analytical approximation and compared with the numerical value of the third iteration, and so on.This approach constitutes one of the tests implemented to validate the code in this study.
To converge the code, we proceeded as follows: since the finite difference approximation method used in the Taylor series development is of second order, we postulate that the error is also of second order at ∆x: where fij is the analytical solution, and Fij is the numerical solution obtained, and L2 = Error.Figure 2 shows the two current wave models used to perform this convergence.That is the analytical solution and the numerical solution obtained.
The waveform generated with the Heidler function [8,9,10] was selected as the analytical solution and compared with the numerical solution of the telegrapher's equation to estimate the convergence error.Figure 3 shows a line of average slope m=1.9 using least squares for linear regression, with differential steps of 3 microseconds.The numerical solution converges to the analytical solution to the second order of approximation.ISSN-E: 2542-3401, ISSN-P: 1316-4821 Universidad, Ciencia y Tecnología, Vol. 28, Núm.122, (pp.7-16) The particular case of a lossless line (R=G=0) leads to a differential equation identical to the string equation and whose algorithm is straightforward to implement numerically.
So that you know, the solution of equations ( 10) and ( 11) can be used to verify the convergence of the code in the limits of R=G=0.To verify the validity of the solution of the telegrapher's equation ( 3), in codes where the equation is solved with R=G=0, the current's numerical solution coincides with that of equation (11).A mesh refinement was performed for N=1000, 10000, 200000, and 1000000, and the code remained stable and convergent for all meshes, only with the centered finite difference scheme.With the forward and backward finite difference schemes, the code is unstable for coarse, medium, and fine meshes.

IV. RESULTS
Another test and validation of the code elaborated in MATLAB was a series of tests performed on the code by varying the parameters.A series of plots are shown for G=0 and G≠0 for different numbers of points to change the size of the spatial differentials and test the previously exposed stability criteria.The fact that the conductance parameter is zero implies non-real solutions since there are losses in the real lines.Meanwhile, when G≠0, the solutions are closer to reality.In the case of Figure 4, for several points N≥1000, the code starts to break the CFL condition, and oscillations are appreciated, representing an atypical behavior for a phenomenon of these characteristics.For a N≤200, the CFL condition (ratio between Δt and Δx) is violated.Once the implemented numerical solution was compared with the analytical solution and the codes were validated, we simulated different models.Figure 6 shows the current calculated from Equation (3) for different values of G.In all cases, the start-up time and the differentials were the same.What differentiates them is the attenuation and the stability of the code.

CONCLUSIONS
It has been possible to create software in Matlab designed to be executed in Windows environments, which has allowed an exhaustive analysis of the stability and convergence in the numerical resolution of the Heaviside equation.The results have shown that this software is a reliable tool since it demonstrates stability and convergence of second order, which is consistent with the finite difference scheme applied in discretizing the Heaviside equation.A prominent feature of the software is its ability to maintain the CFL condition during the evolutionary process, thus guaranteeing reliable and consistent results.
This study is valuable to simulations and numerical solutions of partial differential equations.One of the most outstanding aspects lies in the meticulous validation and rigorous tests to which the developed code has been subjected.Specifically, the stability and convergence tests have yielded highly satisfactory results, supporting the efficiency and reliability of the software in question.
It is essential to note that the convergence criteria, a fundamental part of this analysis, are intrinsic to the nature of the simulated physical system and its actual behavior.These criteria are intricate and are significantly influenced by the experience of the simulator or the research team in charge regarding their specific knowledge of the problem in question.Proper consideration of these criteria not only supports the accuracy of the results obtained but also provides a deeper insight into the fidelity of the simulation to the real system.
Finally, this work has not only resulted in high-performance and robust software for the numerical solution of differential equations but has also provided several valuable lessons on the importance of validation, stability, and convergence testing, as well as the influence of convergence criteria on the quality and reliability of the results.These findings add to the continued growth of the field and open doors for future research and developments in numerical simulations.

Figure 3 .
Figure 3. Result of the average line using least squares for convergence.

Figure 4 .
Figure 4. Numerical solution for different values of the number of mesh points N and the conductance parameter, G=0.

Figure 5
Figure5compares the waves obtained numerically with the code in Matlab (both for G=0 and G≠0) and the analytical solution.It can be noticed the difference with G=0.This is because the discharge channel's dissipative characteristics are not considered.The solutions closest to reality are obtained when G≠0.

Figure 5 .
Figure 5.Comparison between the numerical solutions for the conductance parameter G=0 and G≠0 and the analytical solution (Heidler).This plot was obtained with a run of the program in Matlab..

Figure 6 .
Figure 6.Numerical solution for the discharge current with different values of the conductance parameter G..

Figure 7
Figure7shows a series of plots for different values of the discharge channel resistance, R, both for G=0, to simulate the different resistivities present in the atmosphere by having areas with higher humidity and pollution than others.

Figure 7 .
Figure 7. Numerical solution with different values of channel resistance R and conductance parameter G=0.