Improving Numerical Accuracy in Time-Domain Simulation for Power Electronics Circuits

In time-domain simulations of power system transients, trapezoidal integration with fixed step-size is the most common method due to its accuracy and ease of implementation. Discontinuities occurring within fixed time-step when simulating power electronics circuits, may cause numerical oscillations and errors. Several methods are available in the literature for interpolation and handling of discontinuities. This paper intends to analyze how accuracy is affected by existing techniques for handling discontinuities in time-domain simulations based on the trapezoidal integration method. New algorithms are proposed to improve accuracy.

T HE trapezoidal integration method is a one-step implicit A-stable method commonly employed in the simulation of electromagnetic transients [1], [2]. It has a second order accuracy level [3], [4] the same as 2S-DIRK [5], [6] and Gear 2 [7] and can perform efficiently/accurately for a wide range of transient frequencies [4]. But this method is not perfect.
Firstly, trapezoidal integration causes numerical oscillations at discontinuities. Discontinuities occur due to switching in power electronics circuits or due to nonlinear functions. Several techniques are proposed in the literature to eliminate numerical oscillations. One of them consists of switching from trapezoidal (TR) integration to L-stable Backward Euler (BE) method for two half time-steps after the discontinuity occurrence [4], [8]. This method is implemented in [9]. It follows the initial idea of [10] where the BE method is used only to reinitialize TR integration after locating the discontinuity instant.
Some works have demonstrated the effectiveness of the L-stable 2S-DIRK [5], [6], [11] method for eliminating numerical oscillations. This 2-step method is more complicated to implement than TR and may become less accurate for very stiff systems [4]. TR integration with 2S-DIRK for eliminating numerical oscillations is also tested in [4].
Another method for avoiding numerical oscillations is the chatter removal technique [12]. It consists in performing a half time-step interpolation after crossing a discontinuity [12].
In [13] TR integration is combined with the backward differentiation formula method, called BDF2 [14] to obtain the TR-BDF2 formula.
All of the above approaches may affect numerical accuracy, especially with power-electronics circuits.
Secondly, with fixed time-step simulation, discontinuities occurring between two discrete time points may not be accurately tracked, resulting in additional errors and even in non-characteristic harmonics for some power-electronic circuits [15], [16].
Interpolation techniques are used in [17], [18] to fall back to discontinuity instant and to resynchronize with simulation time-mesh. That is also called double interpolation [12]. This mechanism is assumed to yield better accuracy, compared VOLUME 8, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ to performing fixed step-size numerical integration without applying interpolation. However, sometimes this might not be the case, as it was pointed out in [3] without further analysis.
This paper follows research initiated in [3] to present new details and analysis on accuracy for existing numerical methods with techniques to handle discontinuities. Contrary to [3] the focus of this paper is on TR integration only, since it is commonly used in EMT-type software packages [1]. In addition, this paper presents a more comprehensive analysis on a larger set of existing numerical methods and provides new explanations on accuracy issues while proposing better methods. Also, this paper uses more sophisticated test cases for pin-pointing all potential accuracy issues in switching systems. The contributions are valuable for EMT-type software development and accuracy analysis.
This paper is organized as follows. Section II presents existing techniques for handling discontinuities. Section III contributes new methods and alternatives. In section IV, numerical results based on analysis performed in section II, are presented. Finally, section V shows simulations for practical cases.

II. TECHNIQUES FOR HANDLING DISCONTINUITIES WITH TRAPEZOIDAL INTEGRATION METHOD A. ACCURATE TRACKING OF DISCONTINUITIES
In fixed time-step t computations, when a discontinuity occurs between two consecutive time-points, interpolation can be used to calculate unknown variables at the discontinuity instant. Interpolation can be linear or quadratic. In the case of linear interpolation: where x n , x n+1 and x z are respectively the values of a variable x at time-points t n , t n+1 = t n + t and at the discontinuity instant t z ∈ [t n , t n+1 ]. Quadratic interpolation [3], [19] requires the solution x n−1 at t n−1 and gives: where In [17], [20] and [21] after linear interpolation for solution at t z using (1), trapezoidal integration is performed to reach t z + t. After this point, another linear interpolation is used to resynchronize with time-mesh solution at t n+1 . This mechanism is hereinafter referred to as TR_DI for trapezoidal rule and double interpolation. Resynchronization remains optional and the single interpolation approach is given the acronym TR_SI.
In forced or naturally commutated power electronics circuits, the change of state of one switch may lead to changes of states in one or more switches. Accurate tracking of such discontinuities requires to complete all switch status changes at the same simulation time-point (without advancing in time). This method was initially introduced in EMTP [9] (see also [2]) and named simultaneous switching (SS). In fact, completely wrong answers may result if SS is not applied in some circuit topologies, such as the buck-boost converter.

B. SUPPRESSION OF NUMERICAL OSCILLATIONS
In time-domain simulation, network components are described by differential equations [22] in the form: TR integration can be used to discretize (4): Equation (5) is used to create companion models for all network components and formulate network equations using, for example, modified-augmented-nodal-analysis [9]: where A t is the simulated network Jacobian matrix, x t is the vector of unknowns and b t is the vector of known quantities (includes history terms and known sources). Equation (6) is solved at each time-point t. Various discontinuities may occur due to switching devices in power-electronics circuits or nonlinear functions. It can be shown that although the trapezoidal rule is A-stable, such discontinuities may cause numerical oscillations in (5). Numerical oscillations can be eliminated using various techniques, one of them is presented in [10]. It uses interpolation and a single t/2 projection step BE integration for reinitialization of TR integration.
It is recalled that BE numerical integration is given by This approach [10] can however produce significant errors for large values of t. It is modified in [8] by applying two consecutive t/2 time-step BE integrations without interpolation. This version is used in [9] where the BE steps are applied until all discontinuities are eliminated and followed by TR integration. This method is hereinafter referred to as TR_BE. In [9] discontinuities are automatically detected by monitoring events, such as switch status changes, sources starting and stopping and nonlinear functions changing operating segments. In addition, [9] uses a simultaneous switching algorithm, as also proposed in [21]. Simultaneous switching consists in verifying all switch states before advancing in time.
It is also possible to apply a two-step integration method. The first one is performed with TR and the second one uses an L-stable method. In TR-BDF2 [3], [13], [14], [23],an intermediate stage x γ is calculated to move from x n to x n+1 as: 158 VOLUME 8, 2021 where γ = 2 − √ 2 and η = (1 + √ 2)/2. The internal trapezoidal stage has the size γ t, so t γ = t n + γ t. This idea is similar to the 2S-DIRK method [11] where an intermediate point solution is also used. Such methods are typically more complicated to implement in large scale simulation tools since they require two history term calculations per simulation time-point. This also creates complications for several models, such as rotating machines and transmission lines. Comparative efficiency with TR is a concern. Nevertheless, the target of research in this paper is the TR method.
It is also possible to suppress numerical oscillations by applying a chatter removal technique [12], [17]. This method can be explained using the companion model of an inductance. The voltage across an inductance L at any time-point t n+1 is given by when using trapezoidal integration. When the current i n+1 = 0 due to a discontinuity (switch opening), the voltage becomes v n+1 = −(2L/ t)i n − v n . At the following timepoint i n+2 = 0 and v n+2 = −v n+1 . At any following timepoints the voltage is the opposite of the preceding time-point which constitutes numerical oscillations. The oscillations can be avoided by performing a half time step linear interpolation between time-points t n+2 and at t n+1 . The resulting value (theoretically exactly zero) will stand for the solution at t n+2 and the simulation can continue from there. In practice, some simulation techniques combine double interpolation and chatter removal to handle discontinuities that occur inbetween time-points [12], [17]. It is assumed in this method that perfect equality is reached between two consecutive opposite voltage values, which remains true within negligible numerical roundoff error. We refer to this method as TR_DICR. Note that TR_SI and TR_DI do not remove numerical oscillations.

C. ERROR ANALYSIS
In time-domain simulations with TR integration, the techniques used to handle discontinuities may cause additional errors. This section aims to quantify those errors. The localerror e(t k ) of a given variable x at a given timepoint t k is evaluated by calculating the difference between its numerical value x k and its exact value x(t k ): The exact value x(t k ) may be obtained through a reference simulation or can be analytical. From (11), relative-localerror is defined as: Local and relative-local errors ( (11) and (12)) may not be sufficient to assess accuracy. For a simulation between timepoints t start and t end , the error should be evaluated for every instant t k ∈ [t start , t end ]. To do that, we define the error function for every t k and quantify the overall error by the relative rms error (RRE) as presented in [3]: where x rms is the rms value of the exact solution waveform and n s is the number of simulation points. The local-error of a numerical integration technique after crossing a discontinuity can be evaluated using the basic function [22]:ẋ = λx (14) where Re(λ) < 0. If x 0 is the initial value, the exact solution of (14) is: In the timeline of Fig. 1 the hatted variables are on the time-mesh and the simulation starts fromx 0 att 0 , assuming thatt 0 = 0. A discontinuity occurs at t 01 and α is calculated using (2). According to (14) and (5),x 1 is found in a first step by TR integration:x where . With TR_SI the second step is a linear interpolation to find x 01 usingx 0 andx 1 . According to (1) and (16) The third step is to apply trapezoidal integration to move from x 01 to x 12 : Combining (18), (17) and (16), then using (11) and (15), the e(t 12 ) for TR_SI is given by If we use TR_DI, there is one more step that consists of applying a second interpolation to recalculate the solutionx 1 att 1 , using x 01 and x 12 : From (20), (18), (17) and (16), using (11) and (15), the e(t 1 ) with TR_DI is Now if we use TR_DICR in one of its two forms presented in [12], then from (20), a one-step integration is done to calculate the solution att 2 , that is:x 2 = K 0 K 3 x 0 . Then a half time step interpolation is performed fromx 1 andx 2 . The final solution att 2 becomes: which gives for e(t 2 ) Another approach for performing chatter removal, as explained in [12], consists in performing a half timestep interpolation from x 01 and x 12 in (18), to obtain x 11 . Then a one-step integration is used to calculate x 21 . Finally, linear interpolation is performed from x 11 and x 21 to getx 1 . This approach executes the same steps as TR_DICR, but in a slightly different order. Simulation results not shown here indicate that those two approaches yield similar results, the main objective being avoiding numerical oscillations.
Finally, if we use TR_BE, a discontinuity occurring at t 01 is detected att 1 . The error due to inaccurate detection is the difference betweenx 1 in (16) and the exact solution at t 01 , because both represent instants at which the discontinuity is taken into account. The e(t 1 ) is therefore equal to: Equations (19), (21), (23) and (24) indicate that the instant at which a discontinuity occurs (defined by α) impacts on e. Furthermore, through K 0 , the errors also depend on λ (which characterizes the simulated system) and t. The above equations allow to compare different techniques.

III. NEW METHODS
Finding a perfect numerical method capable to simulate power-electronics circuits at the best accuracy level is a challenging task. Several techniques are developed and used in existing software. They all seek the optimal trade-off between accuracy and computing time. In this section, some new simulation approaches derived from existing methods are proposed. They aim to improve accuracy in simulation tools.

A. TR_BE WITH LINEAR INTERPOLATION (TR_BE_I)
It is possible to improve TR_BE by performing a linear interpolation to find x 01 , before applying the two BE steps. We refer to this method as TR_BE_I. It has not been previously applied in the literature. According to Fig. 1, with TR_BE_I, the simulation moves fromt 0 tot 1 and comes back to t 01 . The result is (17). A first half time-step BE integration is performed to move from t 01 to t 11 . When (7) is applied to (14) x 11 = K 5 x 01 (25) where Then, a second half time-step BE integration is performed to calculate the solution at t 12 : From t 12 , TR integration is used to move the simulation forward without resynchronization to time-mesh. (26) with (17) gives The e(t 12 ) is therefore given by:

B. TR_BE_I WITH TWO TIME-STEP VALUES (TR_BE_I_V)
One way to achieve higher levels of accuracy in simulations is to vary the time-step. The method proposed here is following the work of [24]. We refer to it as TR_BE_I_V. It uses two time-step values: a normal t 1 and a very smaller t 2 . The simulation is normally performed using t 1 . When a discontinuity occurs, linear interpolation is performed to bring the simulation back to the point of discontinuity. At this point the time-step is set to t 2 , two halved time-steps BE integrations are performed, and the time-step is set back to t 1 to continue with TR. This method is similar to TR_BE_I with the difference that the two halved time-step BE integrations are performed using a very small time-step t 2 . The objective is to increase accuracy near the discontinuity point solution.
This method requires to modify A t in (6) for changing the time-step, leading to higher computational burden. Although two versions of A t can be stored and reused.

C. MODIFIED TRAPEZOIDAL WITHOUT INTERPOLATION (MTR_BE_V)
This method allows to directly reach any discontinuity point from the previously calculated time-point by modifying the integration method. Therefore, it accurately tracks discontinuities without interpolation. It is inspired from [25] and [26] where a modified version of the trapezoidal method called weighted-averaged integration was presented as where w is a weight that modifies the numerical method. w = 0 corresponds to BE and w = 1 corresponds to TR. In this formulation, it is possible to substitute t by γ t. Then, by varying γ between 0 and 1, the step-size of the obtained numerical method is varied. If w is adjusted in a way that the condition γ (2 − w) = 1 is always satisfied, then A t in (6) remains unchanged. The new formulation becomes: This new formulation is referred to as MTR_BE_V. It works as follows: the simulation is normally done with γ = 1 (TR integration). At each time-point t n , the code verifies if a discontinuity will occur between t n and t n+1 by checking control signals or using linear extrapolation to make a prediction, as described in [27]. If so, the discontinuity instant t z is calculated and the three following steps are executed: 1) set γ = (t z − t n )/ t and perform one integration step using (30) to advance from t n to t z . 2) set γ = (t n+1 − t z )/ t, account for the discontinuity, then use (30) to advance from t z to t n+1 . 3) set back γ = 1 and perform two halved time-step BE integrations, then continue with trapezoidal integration. This step eliminates numerical oscillations. According to the analysis performed in section II.C, the e(t 2 ) for the proposed MTR_BE_V method, is given by: It is noted that the implementation of MTR_BE_V requires the use of linear extrapolation. This could lead to loss of accuracy as it will be shown later is section IV.B.

D. METHODS WITH SIMULTANEOUS SWITCHING (TR_BE_SS, TR_BE_SS_I)
The simultaneous switching (SS) technique discussed in section II.A. can be combined to TR_BE to obtain the TR_BE_SS method. According to Fig. 1, the simulation moves fromt 0 tot 1 , where the discontinuity is detected. The switch state is changed, and the circuit is solved again att 1 . The algorithm continues to monitor all switches and resolves the circuit until no more changes are detected. The next steps fromt 1 are the same as in TR_BE.
The SS technique can be also included with TR_BE_I to give the method named TR_BE_SS_I. It consists of applying SS att 1 , the same way as in TR_BE_SS and then interpolate to t 01 . The next steps from t 01 are the same as in TR_BE_I. This new method performs both interpolation and SS.
SS is not required for all circuit simulation cases, but when needed, it allows to achieve more accurate results.
The accuracies of above methods will be evaluated next.

A. ACCURACIES OF TECHNIQUES USED TO HANDLE DISCONTINUITIES
In this part, the accuracies of techniques presented to handle discontinuities are evaluated. Equations (19), It is observed from Fig. 2, that without any interpolation technique (TR_BE) to track discontinuities, |e| is linearly correlated with α. TR_SI and TR_DI both use interpolation,  but TR_SI which does not interpolate back to resynchronize with the time-mesh, yields lower errors than TR_DI (uses resynchronization). Therefore, using interpolation for time-mesh resynchronization adds extra errors and should be avoided when unnecessary. It is noted that TR_SI is not capable of removing numerical oscillations and should not be used for practical applications.
With TR_DICR the error does not vary significantly with respect to α. However, it is observed that for higher values of α, TR_BE (no interpolation) performs better than TR_DICR. The above explains the findings of [3] where only comparisons were made without further analysis.
Finally, TR_BE_I and the proposed MTR_BE_V both give lower errors than TR_DICR and TR_BE.
The results presented above are related to local error after discontinuity. In practical simulations, this error may become attenuated in the following steps. It is important to check how accuracy is affected in a complete simulation to verify if the conclusions drawn from Fig. 2 remain valid. To do that, the simple example of Fig. 3, taken from [3] with slight modifications, is simulated with t = 10µs.
The switch SW1 is closed at different instants t z between two consecutive time-points t n = 30ms and t n+1 = t n + t. This means that α defined in (2) varies from 0 to 1 with a step of 0.1: t z = 30ms + α t. The methods TR_BE, TR-BDF2, TR_DICR, MTR_BE_V and TR_BE_I are tested for each value of α. For α = 0 and α = 1, the TR_DICR and TR_BE_I methods account for exact discontinuity instant t z , whereas TR_BE and TR-BDF2 detect it at t n+1 . For each simulation, the relative-rms-error (RRE) in the inductance L1 current is calculated using (11) and (13) for the duration of the transient state from t n+1 = 30ms + t to 40ms. The exact solution in (11), is found using t = 1µs. The results are presented in Fig. 4.  These results validate the theoretical analysis presented in Fig. 2. Indeed, for TR_BE the RRE is again linearly corelated with α. For TR_DICR and TR_BE_I, the RRE is almost constant with varying α. One very interesting observation is the fact that the TR_BE and the TR_DICR graphs intercept each other for a certain value of α. This again explains the counter intuitive result obtained in [3]. When handling a discontinuity, TR_BE may become more accurate than TR_DICR and that depends on α. In a practical simulation, α is unpredictable, so accuracy comparisons amongst some methods may become probabilistic. If SW1 is closed and opened several times, we can claim that TR_DICR that uses interpolation has only about 50% of chance to be more accurate than TR_BE. Also, TR_DICR has 80% more chances to be more accurate than TR-BDF2. TR_BE performs better than TR-BDF2 and the proposed MTR_BE_V performs better than TR_DICR. Furthermore, TR_BE_I is the most accurate technique

B. COMPARISONS BETWEEN INTERPOLATION AND EXTRAPOLATION
Interpolation or extrapolation [27] may be used in simulation algorithms to move back to a discontinuity instant or to calculate the exact moment at which the discontinuity has occurred. We consider 4 consecutive time-points t n−2 , t n−1 , t n and t n+1 in a simulation and assume a discontinuity occurring between t n and t n+1 . To determine the exact instant t z , one can apply linear interpolation from solutions at t n and t n+1 , quadratic interpolation from t n−1 , t n and t n+1 , linear extrapolation from t n−1 and t n or quadratic extrapolation from t n−2 , t n−1 and t n .
This section aims to compare accuracy for the above interpolation and extrapolation techniques. For this purpose, the simple diode rectifier circuit of Fig. 5 is simulated with TR_BE. Each diode is modelled by a 1m resistance in series with an ideal switch model and with a 0.7V DC source when closed, and by a 1M resistance in parallel with an ideal switch when opened.
The circuit is first simulated with t = 100ns (fixed time-step) using TR_BE as reference. Each time a diode state changes, the switching instant t ZR is recorded, no matter what diode it is. t i ZR is the ith value of change. In the second step, the circuit is simulated with TR_BE and with different t values.
For each simulation, each time a diode state changes, linear interpolation, quadratic interpolation, linear extrapolation and quadratic extrapolation are successively used to calculate and record the exact instant. Each obtained value of t Z , named t i Z , is compared to its corresponding value t i ZR from the reference simulation, then the absolute value of the mean relative error is calculated using (32). The results are presented in Fig. 6 where solid lines represent interpolation and dashed lines represent extrapolation. The mean relative error formula is The results indicate that linear interpolation is the most accurate. Differences with quadratic interpolation are negligible. Which means that the more complex quadratic interpolation does not improve accuracy. Also, extrapolation as performed in [27] (linear or quadratic) gives larger errors.
At this level, several techniques to handle discontinuities have been presented and evaluated. For the next section, the most promising amongst them are used to simulate actual power electronics circuits.

A. SIMULATION OF A THYRISTOR CONTROLLED REACTOR (TCR)
The circuit presented in Fig. 7 is taken from [3] and used here for comparison purposes. The thyristors are modeled as the diodes for Fig. 5 in on/off states. The reference simulation is obtained with TR_BE with a 1µs step-size. The circuit is simulated using TR_BE, TR_BE_I, TR_DICR and MTR_BE_V with t = 50µs. The error functions for those methods, obtained using (11) at each time-point are presented in Fig. 8. Each time a thyristor closes, the value of α as defined in (2) is calculated and shown in the graph. In cases of thyristor openings, the current in L1 becomes 0 for all methods and α recording is pointless.
From Fig. 8, it is observed that TR_DICR gives an almost constant error when the thyristors are closed. For TR_BE, the smaller values of α lead to larger errors. These two observations are in accordance with the theory presented in Fig. 2 and with simulation results of Fig. 4. Also, the proposed MTR_BE_V is more accurate than TR_DICR. However other simulation results not shown here, demonstrate that the accuracy of MTR_BE_V can decrease when simulating circuits with diodes. This is because linear extrapolation is used in MTR_BE_V to find discontinuity instants. As stated before and shown in Fig. 6, linear extrapolation degrades accuracy.
Finally, it is noticed that TR_BE_I achieves the best accuracy level and remains the best approach to handle  discontinuities. Also, it is worth to notice that the TR_BE_I_V method (not shown here) gives about the same level of accuracy than TR_BE_I for this circuit. The TCR circuit is a simple circuit. A more complex circuit is simulated in the following section.

B. SIMULATION OF A DC-AC-DC CONVERTER
The circuit of Fig. 9 is made up with 4 IGBTs and 4 diodes. Each IGBT is modelled by 2 diodes and one controlled switch. Therefore, there are 16 switches. Each diode is modelled as in the circuit of Fig. 5. This example demonstrates the previously discussed (section II.V-A) issue in the simulation of electronic converters: the change of state of one or some switches leads to the change of state of some other switches at the same instant. That is why the SS method in is applied below.
The controlled switches of the 4 IGBTs are driven by a periodic step signal generator with 50% width and a frequency of 300 kHz. Five methods are tested: TR_BE, TR_BE_I, TR_BE_I_V, TR_BE_SS and TR_BE_SS_I.
The reference for comparisons is a TR_BE_SS simulation with t = 1 ns. All tested methods are simulated using t = 80 ns. The proposed TR_BE_I_V method uses 80ns and 10ns time-step values. Simulation results for voltage across load R load and absolute error in the same voltage are presented in Fig. 10 and Fig. 11.    Fig. 11 show that TR_BE yields less accurate results for the given time-step. Secondly, we can see that TR_BE_I is more accurate than TR_BE but is less accurate than TR_BE_SS. This reveals that using interpolation alone is not enough to accurately simulate circuits with possible simultaneous commutations. In [21] and [28], it is demonstrated that using a technique to account for instantaneous switching, can reduce the observed losses in power electronics converters. Thirdly, TR_BE_SS_I is the most accurate, compared to TR_BE, TR_BE_I and TR_BE_SS.
The proposed TR_BE_I_V method gives the best accuracy. The very small time-step used after interpolation helps to avoid possible errors due to instantaneous switching situations. It is noted that reducing t improves accuracy for all methods, by shifting them towards the reference solution and in the same order. If the simulation time-step is sufficiently small, all methods deliver similar accuracy. This can be observed in Fig. 12 where the methods are simulated using t = 10 ns. Simulation results not shown here, demonstrate that when the frequency of the control signal generator varies, the conclusions remain the same.

C. EFFICIENCY ASSESSMENT
This part aims to assess simulation efficiency. Different algorithms exhibit different behavior. For instance, the number of times N S that (6) is solved, may vary from one algorithm to another. When simulating the circuit of Fig. 9 for 1 ms, the steady-state error (SSE) and N S are presented in Table 1.
From Table 1, it is seen that the number of system solutions N S is directly corelated with accuracy. More accurate methods exhibit more system solves. TR_BE_SS and TR_BE_SS_I use simultaneous switching and increase the computational burden. Also, since TR_BE_I_V sometimes uses a smaller time-step, the number of simulation points required to complete the simulation increases. Therefore, the computation burden increases.
For the 4 methods compared above, we can normalize criteria values by dividing them by the maximum of the 4 values in each row (see Table 2).
The best tradeoff between accuracy and computational burden is obtained for a method with the smallest SSE and N S as expressed by the last row in Table 2. The smallest value which is 0.434 is achieved with TR_BE_I_V, which clearly demonstrates that this method is the most efficient simulation technique for the studied example. It gives the best tradeoff between accuracy and computational burden. Also, it is seen from the last row of Table 2 that TR_BE_SS gives a better tradeoff than TR_BE_I. In this case, usage of simultaneous switching technique is more efficient than using only interpolation. Finally, TR_BE_SS_I gives the best tradeoff, compared to TR_BE_I and TR_BE_SS.

VI. CONCLUSION
In this paper, we presented new analysis on how interpolation affects accuracy with trapezoidal integration. We demonstrated that interpolation may lead to different accuracy levels, depending on when exactly the discontinuity occurs and depending on the technique used to handle it. We demonstrated the limits of interpolation by showing and explaining why: a) a method with no interpolation can sometimes be more accurate; b) interpolation alone cannot be sufficiently accurate when simulating simultaneously switching circuits.
Usage of Backward Euler method combined with trapezoidal integration, interpolation and variable time-step (TR_BE_I_V), provides the most accurate results. However, due to related complexity in software programming with TR_BE_I_V, it can be concluded that the fixed time-step Backward Euler method combined with trapezoidal integration, interpolation and simultaneous switching (TR_BE_SS_I) is the best tradeoff between complexity and accuracy.