ISSN-E: 2542-3401, ISSN-P: 1316-4821

Universidad, Ciencia y Tecnología,

Vol. 28, Núm. 122, (pp. 7-16)

Suárez et al. Convergence and stability criteria for numerical solutions of partial differential equations in science and engineering

https://doi.org/10.47460/uct.v28i122.761

Convergence and stability criteria for numerical

solutions of partial differential equations in

science and engineering

Franyelit Suárez-Carreño

http://orcid.org/0000-0002-8763-5513

franyelit.suarez@udla.edu.ec

Universidad de las Américas

Facultad de Ingeniería y Ciencias Aplicadas

Carrera de Ingeniería Industrial

Quito, Ecuador

Luis Rosales-Romero

https://orcid.org/0000-0002-7787-9178

luis.rosales2@gmail.com

Universidad Nacional Experimental Politécnica

Antonio José de Sucre

Vice Rectorado Puerto Ordaz,

Doctorado en Ciencias de la Ingeniería

Puerto Ordaz, Venezuela

Recibido (14/08/2023), Aceptado (23/11/2023)

Abstract. - 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.

Keywords: simulation, convergence, stability, atmospheric discharges.

Criterios de convergencia y estabilidad para soluciones numéricas de ecuaciones

diferenciales parciales en ciencia e ingeniería

Resumen: En este trabajo se presenta una exploración de los aspectos informáticos de la simulación y el

modelado aplicados a la resolución de la ecuación de Heaviside, considerando la relevancia de la estabilidad y

convergencia de las soluciones. Para ello se implementó un esquema de diferencias finitas de segundo orden

como enfoque principal para el estudio de las descargas atmosféricas. El lenguaje de programación utilizado

fue Matlab, lo que facilitó el cálculo de las corrientes inducidas en el escenario de estudio. Para la

implementación numérica se consideraron los enfoques de diferencia centrada, diferencia hacia adelante y

diferencias finitas hacia atrás. La validación del sistema se hicieron pruebas que demuestran la efectividad del

diseño y la convergencia al segundo orden, con el enfoque de diferencia centrada.

Palabras clave: simulación, convergencia, estabilidad, descargas atmosféricas.

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

ISSN-E: 2542-3401, ISSN-P: 1316-4821

Universidad, Ciencia y Tecnología,

Vol. 28, Núm. 122, (pp. 7-16)

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:

Suárez et al. Convergence and stability criteria for numerical solutions of partial differential equations in science and engineering

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.

ISSN-E: 2542-3401, ISSN-P: 1316-4821

Universidad, Ciencia y Tecnología,

Vol. 28, Núm. 122, (pp. 7-16)

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

Suárez et al. Convergence and stability criteria for numerical solutions of partial differential equations in science and engineering

A. The Algorithm

The algorithm for finding the numerical solution is described in Figure 1.

ISSN-E: 2542-3401, ISSN-P: 1316-4821

Universidad, Ciencia y Tecnología,

Vol. 28, Núm. 122, (pp. 7-16)

Figure 1. The algorithm for finding the numerical solution.

Source: Authors.

Equation (4) is implemented numerically in MATLAB. In the case of the telegrapher's equation, the input

parameters R, L, C, and G are needed, which are assumed to be known in the discharge channel. The values

used in this case are shown in Table 1.

Table 1. Electrical parameters of the discharge channel of an ATDI. These parameters were

varied to validate the code, and the impact on the solutions was not representative.

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.

Suárez et al. Convergence and stability criteria for numerical solutions of partial differential equations in science and engineering

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:

ISSN-E: 2542-3401, ISSN-P: 1316-4821

Universidad, Ciencia y Tecnología,

Vol. 28, Núm. 122, (pp. 7-16)

By applying logarithm and its properties to both members of the equation (7):

Equation (8) has a line of slope equal to two. If from the graph Vs the error between a known analytical

solution and the numerical solution found for equation (4) is obtained. A slope approximately equal to two

implies that the numerical solution converges to the analytical solution. In this case, a general expression for

the convergence calculation given by the L2 norm was used for convergence [3,4,11]

Suárez et al. Convergence and stability criteria for numerical solutions of partial differential equations in science and engineering

ISSN-E: 2542-3401, ISSN-P: 1316-4821

Universidad, Ciencia y Tecnología,

Vol. 28, Núm. 122, (pp. 7-16)

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.

Figure 2. Numerical and analytical discharge current waveforms.

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

Suárez et al. Convergence and stability criteria for numerical solutions of partial differential equations in science and engineering

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.

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

Suárez et al. Convergence and stability criteria for numerical solutions of partial differential equations in science and engineering

ISSN-E: 2542-3401, ISSN-P: 1316-4821

Universidad, Ciencia y Tecnología,

Vol. 28, Núm. 122, (pp. 7-16)

Figure 5 compares 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. 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..

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.

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

Suárez et al. Convergence and stability criteria for numerical solutions of partial differential equations in science and engineering

ISSN-E: 2542-3401, ISSN-P: 1316-4821

Universidad, Ciencia y Tecnología,

Vol. 28, Núm. 122, (pp. 7-16)

Figure 7 shows 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. Numerical solution with different values of channel resistance R and conductance parameter G=0.

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.

Suárez et al. Convergence and stability criteria for numerical solutions of partial differential equations in science and engineering

ISSN-E: 2542-3401, ISSN-P: 1316-4821

Universidad, Ciencia y Tecnología,

Vol. 28, Núm. 122, (pp. 7-16)

REFERENCES

[1] R. Álvarez and L. Rosales “Simulación de descargas atmosféricas y su efecto en líneas eléctricas de

potencia”, UCT, V.16 Nro 64 2011.

[2] P. Rodríguez, J. Toledo, E. Contreras, and L. Rosales “Simulación de descargas atmosféricas mediante la

ecuación de onda viajera”, UCT V. 13, Nro 53 2009.

[3] L. Rosales, W. Barreto, C. Peralta, and B. Rodríguez “Nonadiabatic charged evolution in the postquasistatic

approximation” Phys. Rev. D, 40, 24 2010.

[4] Thomas, J., “Numerical Partial Differential Equations-Finite Difference Methods,” USA, Springer-Verlag New

York, Inc., 1995, pp. 205-258.

[5] Hwang, C., “Numerical Modeling of Lightning Based on the Traveling Wave Equations”, IEEE Transactions on

Industry Applications, Vol. 33, Nº. 2, 1997, pp. 1520-1523

[6] McCann, G. et al., “Lightning Phenomena, Chap. 12, Electric Transmission and Distribution Reference Book”,

Pittsburg, Westinghouse Electric & Manufacturing Co., 1950.

[7] Torres, H., “El Rayo, Mitos, Leyendas, Ciencia y Tecnología”, 2daEdición, Bogotá, UNIBIBLOS, 2002.

[8] Heidler, F. et al., “Calculation of Lightning Current Parameters”, IEEE Transactions on Power Delivery, Vol. 14,

N º. 2, 1999, pp. 309-404.

[9] Herrera, J., “Nuevas Aproximaciones en el Cálculo de Tensiones Inducidas por Descargas Eléctricas

Atmosféricas”, Programa de Doctorado en Ingeniería Eléctrica, Facultad de Ingeniería, Departamento de

Ingeniería Eléctrica y Electrónica, Universidad Nacional de Colombia, Colombia, 2006.

[10] Agoris, D., Charalambakos, V., Pyrgloti, E., and Grzybowski, S. (2004). A computational approach to the

study of Franklin rod height impact on striking distance using a stochastic model. J. Electrostatics. 60: 175-181.

[11] Cole, P., Krehbiel, T.J., Henebry, and G.M. (2016). Web-Enabled Landsat Data Time Series for Monitoring

Urban Heat Island Impacts on Land Surface Phenology. IEEE Journal Selected Topics in Applied Earth

Observations and Remote Sensing, 9: 2043-2050.

[12] Rakov, V., and Uman, M. (2003). Lightning: Physics and effects” Cambridge University Press. Cambridge,

UK.

[13] Torres, H. (2015). El rayo en el trópico. Colección apuntes maestros, Ed. UN, Bogotá, Colombia.

Suárez et al. Convergence and stability criteria for numerical solutions of partial differential equations in science and engineering