Review Article  Open Access
Jorge F. Oliveira, José C. Pedro, "Radio Frequency Numerical Simulation Techniques Based on Multirate RungeKutta Schemes", Journal of Applied Mathematics, vol. 2012, Article ID 528045, 21 pages, 2012. https://doi.org/10.1155/2012/528045
Radio Frequency Numerical Simulation Techniques Based on Multirate RungeKutta Schemes
Abstract
Electronic circuit simulation, especially for radio frequency (RF) and microwave telecommunications, is being challenged by increasingly complex applications presenting signals of very different nature and evolving on widely separated time scales. In this paper, we will briefly review some recently developed ways to address these challenges, by describing some advanced numerical simulation techniques based on multirate RungeKutta schemes, which operate in the onedimensional time and also within multidimensional frameworks.
1. Introduction
In recent years, radio frequency (RF) and microwave electronic circuit simulation has been conducted to an increasingly demanding scenario of heterogeneous broadband and strongly nonlinear wireless communication circuits, presenting a wide variety of slowly varying and fast changing state variables (node voltages and branch currents). In such circuits, the baseband analog blocks, the digital blocks, and the RF blocks may be all intricately mixed. Classical simulation tools are not capable of handling this kind of circuits in an efficient way because they do not perform any distinction between nodes or blocks within the circuit, considering their time constants and/or excitation regimes. Thus, all the blocks in the circuits are treated in the same way, which means that the same numerical algorithm is required to simultaneously compute the response of the digital blocks, the baseband analog blocks, and the RF blocks. Taking into account that signals in different blocks (different parts of the circuit) have completely different formats and evolve on widely separated time scales (which may differ from three, or more, orders of magnitude), it is easy to conclude that the application of the same numerical scheme to all the blocks will result in high inefficiency.
To cope with this scenario, some innovative timedomain techniques have been proposed in the literature in the last few years. In order to benefit from the different rates of variation of slowly varying and fastvarying currents and voltages (circuits’ state variables) in different parts of the circuit, several advanced numerical techniques based on multirate RungeKutta schemes have been proposed to operate within onedimensional and multivariate frameworks. Such techniques include diverse circuit partitioning strategies, which allow the simulator to automatically split the network into subcircuits according to the different time rates of change of their state variables. In order to provide a general overview on these simulation techniques, this paper is organized as follows.
After this brief introduction, Section 2 provides some general background material on standard timedomain circuit simulation techniques commonly used for computing the numerical solution of ordinary unirate problems. Then, Section 3 is essentially devoted to the presentation of some mathematical details of the numerical simulation algorithms based on multirate RungeKutta schemes. Finally, Section 4 concludes this paper.
2. TimeDomain Circuit Simulation Fundamentals
2.1. Mathematical Model of an Electronic Circuit
The behavior of an electronic circuit can be described with a system of equations involving voltages, currents, charges, and fluxes. This system of equations can be constructed from a circuit description, using, for example, nodal analysis, which involves applying the Kirchhoff current law to each node in the circuit, and applying the constitutive or branch equations to each circuit element. In general, systems generated this way have the form where and stand for the excitation (independent voltage and current sources) and state variable (node voltages and branch currents) vectors, respectively. represents memoryless linear or nonlinear elements, as resistors, nonlinear voltagecontrolled current sources, and so forth, while models dynamic linear or nonlinear elements, as capacitors (represented as linear or nonlinear voltagedependent electric charges), or inductors (represented as linear or nonlinear currentdependent magnetic fluxes).
The system of (2.1) is, in general, a differential algebraic equations’ (DAE) system. For achieving an intuitive explanation of the mathematical formulation of (2.1), let us consider the basic illustrative example depicted in Figure 1. This circuit is composed of a current source connected to a linear inductance and two nonlinear circuit elements commonly used when modeling semiconductor devices (a nonlinear capacitance and a nonlinear voltagedependent current source). These nonlinearities are assumed as quasistatic [1, 2] and thus are described by algebraic constitutive relations of voltagedependent charge and voltagedependent current. A nodal analysis of this circuit leads to the following system of equations in the capacitor voltage and the inductor current : This system can be seen as a particular case (in ) of the DAE system of (2.1), in which the excitation vector, , and the vector of state variables, , are given by
The DAE system of (2.1) may obviously be written in other forms. For instance, if we apply the chain differentiation rule to the dynamic term of its left hand side, we can obtain or in which is usually known as the mass matrix. If this matrix is nonsingular, then the DAE system of (2.4) may degenerate into the following ordinary differential equations’ (ODE) system: which can be rewritten in the classical form commonly used in the mathematical literature. When is singular, the DAE system of (2.4) will not degenerate into an ODE system, but it is often possible to express it as a set of algebraic equations combined with a set of differential equations of the form (2.7).
From the above, we conclude that, in many cases, electronic circuits may be described by ODE systems instead of DAE systems. For example, if we return to the simple nonlinear dynamic circuit of Figure 1 and rewrite (2.2) as or, in its vectormatrix form, as we can easily see that if , then the mass matrix is nonsingular, and the circuit may be described by an ODE system expressed in the classical form of (2.7), which, in this case, will simply result in
2.2. Transient Analysis: TimeStep Integration
Obtaining the solution of (2.1) over a specified time interval with a specific initial condition is what is usually known as an initial value problem, and computing such solution is frequently referred to as transient analysis. The most natural way to evaluate is to numerically timestep integrate (2.1) directly in time domain. One possible way to do that consists in simply converting the differential equations into difference equations, in which the time derivatives are approximated by appropriate incremental ratios. With this strategy, the nonlinear differential equations’ system of (2.1) is converted into a purely nonlinear algebraic system.
The above formulation derives directly from the intuitive idea that derivatives can be approximated, and thus simply replaced, by finitedifference schemes. Although this technique can be used to compute the transient response of a generic electronic circuit described by (2.1), there is as alternative strategy which is more often employed to find the solution of initial value problems. Such strategy consists in using initial value solvers, as linear multistep methods [2–5] or RungeKutta methods [3–5] (the most popular timestep integrators). Both of these classes of methods can provide a wide variety of explicit and implicit numerical schemes, with very distinct properties in terms of order (accuracy) and numerical stability. However, since a substantial part of the research work reviewed in this paper is based on modern multirate RungeKutta schemes, we will restrict our presentation to only these numerical techniques. Linear multistep methods will not be addressed.
2.3. RungeKutta Methods
In view of the fact that the wellestablished theory of numerical integration is oriented toward the solution of standard ODEs, we will now consider the form of (2.7) for the mathematical description of a circuit’s operation. So, let us consider a generic initial value problem with state variables, expressed in its classical form by the system of (2.7) and the initial condition , that is,
Definition 2.1 (RungeKutta (RK) method). A standard sstage RK method expressed by its Butcher tableau [5]
for obtaining the numerical solution of (2.11) at the time instant is defined as [3, 5]
where
The algorithm defined by (2.13), (2.14) allows the numerical solution at any generic time instant to be evaluated from its previous calculated value . If we have for , then each of the in (2.14) is explicitly given in terms of the previously computed , , and the method is then an explicit RungeKutta method. If this is not the case, then the method is implicit and, in general, it is necessary to solve a nonlinear system of algebraic equations to simultaneously compute all the . In general, any iterative technique (e.g., fixedpoint iteration or NewtonRaphson iteration) may be used to solve the nonlinear system of (2.14).
RungeKutta methods are universally utilized for timestep integrating initial value problems and differ from linear multistep methods in several aspects. Since they present a genuine onestep format, one of their main advantages is that there is no difficulty in changing the steplength in a dynamic timestep integration process (in opposition to multistep methods where considerable difficulties may be encountered when we want to change steplength [5]). It must be noted that small step sizes can provide a good accuracy in the simulation results but may conduct to large computation times. On the contrary, large step sizes will reduce the computation time but will definitely conduct to poorer accuracy. A good compromise between accuracy and simulation time is achieved when is dynamically selected according to the solution’s rate of change. The automatic step size control is based on the estimation of local errors, for which diverse techniques can be used, such as extrapolation techniques or embedded RK formulas. All theoretical and technical details of implementation of these techniques, as well as many other aspects of the RK methods, such as consistency, convergence, order conditions, numerical stability, and so forth, can be seen, for example, in [3] or [5].
2.4. Conventional Transient Circuit Simulation: SPICE
In the above, we have seen that the most natural way of simulating an electronic circuit is to numerically timestep integrate the ordinary differential system describing its operation. So, it should be of no surprise that this straightforward technique was used in the first digital computer programs of circuit analysis and is still nowadays the most widely used numerical method for that purpose. It is present in all SPICE (which means simulation program with integrated circuit emphasis) [6] or SPICElike computer programs.
SPICE was initially developed at the Electronics Research Laboratory of the University of California, Berkeley, in the early 1970s. The real popularity of SPICE started with SPICE2 [6] in 1975, which was a muchimproved program than its original version (SPICE1), containing several analyses (AC analysis, DC analysis, DC transfer curve analysis, transient analysis, etc.) and device models needed to design integrated circuits of that time. Other versions of SPICE have been developed along the years, and today many commercial simulators are based on SPICE. However, its application to RF circuits may cause some problems resulting from the specific behavior of RF systems. To understand this, we must recall that RF signals are typically narrowband. This means that a data signal with a relatively low bandwidth is transmitted at a very high carrier frequency. To simulate a sufficient portion of the data signal, for example, to estimate bit error rate in modern wireless systems, a large number of carrier periods must be timestep integrated, and thus a huge number of time samples is required (a large consumption of memory and computational time).
2.5. SteadyState Analysis
Although some simulation tools focus on transient analysis (SPICElike simulation), the steadystate behavior of the circuits is of primary interest to RF and microwave designers. The main reason for this is that wireless systems are expected to handle a sinusoidal RF carrier modulated by a slowly varying envelope, so that most aspects of system performance are easier to characterize and verify in the carrier’s steady state. For instance, harmonic or intermodulation distortion, noise, power, and transfer characteristics such as gain, or impedance are examples of quantities that are best defined in the frequency domain and thus measured when the circuit is in steady state.
Timestep integration engines, as the ones mentioned above, and which are tailored for finding the circuit’s transient response, are not adequate for computing the steadystate. As discussed above, timestep integration is a numerical technique intended to obtain the solution of an initial value problem, as it evaluates for a set of successive time instants (time steps) given an initial condition . However, if the objective is the determination of the steady state, there is no alternative other than to pass through the lengthy process of integrating all transients, and expecting them to vanish. In circuits presenting extremely different time constants, or high resonances, as is typically the case of narrowband RF and microwave circuits, timestep integration can be very inefficient. Indeed, in such cases, frequencies in steadystate response are much higher than the rate at which the circuit approaches this regime, or, in other words, the ratio between the state variables’ highest and the lowest frequency components is very large. Thus, the number of discretization time steps used by the numerical integration scheme will be enormous, because the time interval over which the differential equations must be numerically integrated is set by the lowest frequency, or by how long the circuit spends in achieving the steady state, while the size of the time steps is determined by the highest frequency component.
If the steadystate response of a circuit consists of generic waveforms presenting a common period, then the circuit is said to be in periodic steadystate. Computing the periodic steady state response of an electronic circuit involves finding the initial condition, , for the differential system that describes the circuit’s operation, such that the solution at the end of one period matches the initial condition, that is, , where is the period. Problems of this form, that is, those of finding the solution to a system of ordinary differential equations that satisfies constraints at two or more distinct points in time, are referred to as boundary value problems. In this particular case, we have a periodic boundary value problem that can be formulated as or, in the classical ODE form, as where the condition is known as the periodic boundary condition.
Solving (2.16) involves computing a numerical solution that simultaneously satisfies the differential system and the twopoint periodic boundary condition. Certainly, there is no shortage of mathematical literature describing methods for solving boundary value problems. However, a particular technique has been found especially useful for electronic circuit problems. This technique is known as the shootingNewton method [7].
2.6. The ShootingNewton Method
The timedomain method most commonly used for numerically evaluating the periodic steadystate solution of an electronic circuit is the shooting method. Shooting solves boundary value problems by computing the solution to a succession of initial value problems with progressively improved guesses of an initial condition, which ultimately results in the steady state. In a circuit’s steadystate simulation, shooting begins by simulating the circuit for one period using some guessed initial condition (generally determined from a previous DC analysis). Then, the computed solution at the end of the period is checked, and if it does not agree with the initial condition, the initial condition is wisely modified. The circuit is then resimulated with the adjusted initial condition, and this process is repeated until the solution after one period matches the initial condition.
In order to provide some mathematical details on the implementation of the shooting method, let us consider (2.16). Now suppose that we want to numerically timestep integrate the differential system in (2.16) with an initial value solver. As stated above, timestep integration is tailored for transient analysis but is inadequate for computing steadystate responses. The problem comes from the fact that we do not know a priori which initial condition must be considered that will lead to the steadystate solution in the period , that is, that will satisfy the periodic boundary condition . So, we are trying to solve a boundary value problem with an initial value solution technique. One possible way to convert the initial value solution procedure into a boundary value problem solver consists of guessing the initial estimate of , or shooting for , timestep integrating the differential system from until , comparing the resulting with , and then wisely updates the initial estimate. So, shooting is an iterative solver that uses an initial value technique to solve a boundary value problem. In the end, it relies on finding the solution of Let us now define and rewrite (2.17) as where is the statetransition function [7, 8]. An easy way to solve (2.18) consists in using the fixedpoint iteration solver, which, in this case, would simply result in However, shooting with the fixedpoint iteration technique is obviously equivalent to integrating the original differential system from until all transients decay. So, it is generally a useless technique because the convergence to the periodic steadystate solution may be extremely slow. A wellknown tactic to accelerate the route to steady state, that is, to accelerate the convergence of (2.18) to its solution, is the socalled shootingNewton technique. As happens with any other shooting technique, shooting Newton is based on guessing initial conditions. However, it can take advantage of the fact that, although electronic circuits can be strongly nonlinear, their statetransition functions are usually only mildly nonlinear. This means that slight perturbations on the initial condition (starting state) produce almost proportional perturbations in the subsequent time states. Taking this into account, it is easy to conclude that (2.18) can be iteratively solved in an efficient way with the Newton’s method, which in this case will lead us to where is the identity matrix. The only entity of (2.20) that is difficult to compute is the Jacobian of the statetransition function (usually referred to as the sensitivity matrix). In order to compute this matrix, we must take into consideration the chain differentiation rule. In fact, since is nothing but the numerical value , with being the total number of time steps in the interval , which depends on the previous value , which, itself, depends on , and so forth, the sensitivity matrix can be given by
It is easy to see that all the matrices in (2.21) can be individually computed along the timestep integration process. For concreteness, let us suppose that a standard RungeKutta method (2.13), (2.14) is being used to perform timestep integration on the consecutive iterations of the shooting method. If, for example, we want to evaluate the first matrix, then we have Now, if we rewrite (2.14) as we obtain with
Although solving (2.20) and computing the sensitivity matrix may involve some extracomputational cost, shooting Newton converges to the steadystate solution much faster than the normal timestep integration procedure (shooting with fixedpoint iteration). This is the reason why it is the timedomain steadystate engine most widely used in circuit simulation.
3. Numerical Simulation Algorithms Based on Multirate RungeKutta Schemes
3.1. TimeStep Integration with Different Step Sizes
As stated above, dynamic behavior of some electronic circuits involves signals with widely separated rates of variation. Analog circuits presenting extremely different time constants, coupled systems of analog and digital networks, or combined technologies of RF and baseband analog (or even digital) blocks in the same circuit are typical examples where the corresponding state variables may evolve according to very distinct time scales. In such cases, node voltages or branch currents presenting slow (latent) and very fast (active) time evolution rates may coexist in the same problem. This phenomenon, in which some of the state variables are varying very slowly (or even being practically constant) within a specific time interval, while other variables exhibit fast variations in that interval, is frequently referred to as timedomain latency [9–15]. Another common situation where timedomain latency can be found refers to purely digital circuits. For example, in largescale integration (LSI) or very largescale integration (VLSI) applications, usually, only a small part of the circuits (a small number of state variables) is active in a certain time interval, whereas the majority is latent (a large number of state variables may possibly remain practically constant within that interval). In summary, many kinds of electronic circuits may exhibit timedomain latency.
As described in Section 2, timestep integration is a conventional technique that is used by SPICElike computer programs for simulating electronic circuits. However, when integrating differential systems whose components (state variables) evolve according to different time rates, one would like to use numerical schemes that do not expend unnecessary work on slowly changing components. In fact, in such cases, traditional timestep integrators, like standard RungeKutta or linear multistep methods—which use the same step size for all system’s components—become inefficient. To cope with this, some modern multirate RungeKutta (MRK) schemes have been proposed in the literature in the recent years [11–14]. These powerful schemes split the differential system of (2.11) into coupled active and latent subsystems with where is the active (fastvarying) statevariable components’ vector, and is the latent (slowly varying) statevariable components’ vector. The active components are then integrated with a small step size (microstep), while the latent components are integrated with a much larger step size (macrostep). The number of microsteps within each macrostep is an integer that will be denoted by , that is to say, . This is illustrated in Figure 2.
(a)
(b)
In summary, the vector of active state variables is calculated at each of the time instants , defined in the fine grid , whereas the vector of latent state variables is evaluated only in the coarse time instant . A general definition of a multirate RungeKutta timestep integrator is presented in the following.
Definition 3.1 (multirate RungeKutta (MRK) method). Consider two RungeKutta methods of and stages, that can be, but do not necessarily have to be, the same, expressed by their Butcher tableaus and : The resulting multirate RungeKutta method for obtaining the numerical solution of the partitioned system of (3.1), using a microstep for the active components and a macrostep for the latent components, is defined as follows [11, 12]: (i)the active (fastvarying) components are given by where with (ii)the latent (slowly varying) components are given by where with
As can be seen, the coupling between the active and latent subsystems is performed by the intermediate stage values and . There are several strategies for computing these values, as the ones suggested in [9, 10], which are based on interpolation and/or extrapolation techniques, or the ones more recently proposed in [11, 12] or [13], which are based on coupling coefficients.
It must be noted that, in general, the partition into fast and slow subsystems, as well as the number of microsteps within a macrostep, may vary throughout the integration process. Technical details of MRK code implementation, such as activelatent partitioning strategies, step size control tools, or even stiffness detection stratagems, are described in detail, for example, in [11] or [15].
3.2. Multitime Formulation
Signals handled by wireless communication systems can be usually described by a highfrequency RF carrier modulated by some kind of baseband information signals, as amplitude and/or phase signals. The general form of an amplitude and phasemodulated signal can be expressed as where and are, respectively, the slowly varying amplitude and phase baseband signals (the complex envelope in a constellation diagram), modulating the fastvarying carrier.
Circuits driven by envelope (amplitude and/or phase) modulated signals, or presenting themselves state variables of this type, are common in RF and microwave applications. Since the baseband signals have a spectral content of much lower frequency than the carrier, that is, because they are typically slowly varying signals, while the carrier is a fastvarying entity, simulating nonlinear circuits—in a computationally efficient way—containing this kind of signals is often a very challenging task. Because the aperiodic nature of the envelope modulation signals obviates the use of any steadystate technique, one might think that conventional timestep integration (SPICElike transient simulation) would be the natural method for simulating such circuits, as long as the time scale of the signals’ slowly varying components is comparable to the larger time constants involved. However, the large time constants of the bias networks determine long transient regimes and, as a result, the obligation of simulating a large number of carrier periods. In addition, computing the RF carrier oscillations long enough to obtain information about its envelope properties is, itself, a colossal task. Timestep integration is thus inadequate for simulating this kind of problems.
In this section, we will introduce a powerful strategy for analyzing nonlinear circuits handling multirate signals, that is, signals containing two, or more, entities running on widely separated time scales. This strategy is suitable to deal with the abovedescribed amplitude and/or phasemodulated signals, as also with any other kind of multirate signals. It uses multiple time variables to describe multirate behavior, and it is based on the fact that multirate signals can be represented much more efficiently if they are defined as functions of two or more time variables, that is, if they are defined as multivariate functions [17, 18]. As we will see, with this multivariate formulation, circuits will be no longer described by ordinary differential systems in the onedimensional time t, but, instead, by partial differential systems.
The advantages of using multivariate representations are easily illustrated by considering a bidimensional problem. Thus, let us consider, for example, an amplitude modulated RF carrier of the form where is an envelope, slowly varying in time, while is a fastvarying RF carrier. As explained, simulating a circuit with this kind of stimulus, and using conventional timestep integration schemes, tends to be highly inefficient because it requires time steps closely spaced in time (for representing each fast RF cycle accurately) and during a very long time window determined by the envelope. So, let us now consider a bivariate representation for , defined as where is the slow envelope time scale, and is the fast carrier time scale. In this particular case, is a periodic function with respect to but not to , that is, Figures 3 and 4 depict the univariate and bivariate forms defined in (3.11) and (3.12), respectively, for a time interval and a rectangular region . An envelope and a carrier frequency of , were considered in this basic illustrative example.
By comparing Figures 3 and 4, it can be seen that does not have as many undulations as , allowing thus a more compact representation with fewer samples. Furthermore, due to the periodicity of in , we know that its plot repeats over the rest of this time axis. Thus, the bivariate form plotted in Figure 4 contains all the information necessary to recover the original univariate form depicted in Figure 3. In order to get a realistic idea of the savings that can be achieved by using the bivariate formulation, let us suppose that 20 points were sufficient to represent the slowly varying envelope signal accurately in the interval, and that 20 points were also used per each carrier cycle. In such case, the total number of samples required to represent the bivariate form will be . On the other hand, since there are 50 carrier cycles in the interval, the number of samples required to represent the univariate form will be . Now, if we consider a realistic RF scenario, with a much higher value for the carrier frequency (e.g., ), then the number of points necessary to represent will considerably increase to , while the number of points necessary to represent will remain exactly the same .
Let us now consider a general nonlinear RF circuit described by the differential algebraic equations’ (DAE) system of (2.1), and let us suppose that this circuit is driven by an envelopemodulated signal of the form (3.11). Taking the above considerations into account, we will adopt the following procedure: for the slowly varying parts (envelope time scale) of the expressions of and , is replaced by ; for the fastvarying parts (RF carrier time scale), is replaced by . This results in bivariate representations for the excitation , and the solution and the application of this bivariate strategy to the system of (2.1) converts it into the following multirate partial differential algebraic equations’ (MPDAE) system [17, 18]:
The mathematical relation between (2.1) and (3.14) establishes that if and satisfy (3.14), then the univariate forms and satisfy (2.1) [17]. Therefore, univariate solutions of (2.1) are available on diagonal lines , along the bivariate solutions of (3.14), that is, may be retrieved from its bivariate form , by simply setting . Consequently, if we want to obtain the univariate solution in the generic interval, due to the periodicity of the problem in the dimension, we will have on the rectangular domain , where represents the remainder of division of by . Equation (3.15) defines the sawtooth path illustrated in Figure 5.
The multivariate formulation can be employed to solve many types of multirate problems. It can be adopted to deal with either weakly or strongly nonlinear regimes, as well as with a large class of multirate signals. The main advantage of this MPDAE approach is that it can result in significant improvements in simulation speed when compared to DAEbased alternatives [16–22].
3.3. Envelope following Methods Using RK Schemes
In order to compute the bivariate solutions, some initial/boundary conditions must be added to (3.14). These are determined by the existence, or not, of periodicity in the signal components. As stated above, a typical case of significant practical interest is the envelopemodulated regime, which leads to an initialboundary value problem. Indeed, envelopemodulated responses to excitations of the form of (3.10) correspond to a combination of initial and periodic boundary conditions in (3.14). This means that the bivariate forms of these solutions can be obtained by numerically solving the following initialboundary value problem [17]: on the rectangle . is a given initialcondition function defined on , satisfying [17], and the periodic boundary condition is due to the periodicity of the problem in the fast carrier time scale.
The envelope transient over shooting [16, 17, 19, 23, 24] is an efficient timedomain method that can be used to obtain the numerical solution of circuits described by initialboundary value problems of the form (3.16). This method is a particular implementation of a general technique that is often referred to as envelope following [25] and consists in replacing the derivatives of the slow aperiodic time scale by finitedifference approximations (e.g., the Backward Euler rule), to then obtain a set of successive boundary value problems with periodic boundary conditions where , and is the period in the periodic fast time scale . This means that once is known, the solution on the next slow time instant is obtained by solving (3.17). Thus, for obtaining the whole solution in the entire domain , a total of boundary value problems have to be solved, with being the number of steps in . With the envelope transient over shooting technique, each of the periodic boundary value problems of (3.17) is solved using the shooting method described in Section 2.6, where RK methods can be used to perform the successive timestep integrations in the consecutive shooting iterations.
3.4. Envelope following Methods Using MRK Schemes
Recently, a very powerful computeraided design tool especially conceived for the efficient timedomain simulation of highly heterogeneous nonlinear RF circuits has been proposed [16, 19]. This technique is based on an ingenious modification of the abovedescribed envelope transient over shooting technique. It splits the circuits into two distinct parts, according to the time rates of change of their state variables, and, instead if using classical timestep integrators, it uses modern multirate RungeKutta (MRK) schemes to perform timestep integration with different step sizes in each of the consecutive shooting iterations needed to solve (3.17). Indeed, the system of (3.17) is partitioned into the following active and latent subsystems: where is the active (fastvarying) state variable components’ vector at the slow time instant , and is the latent (slowly varying) state variable components’ vector at the same slow time instant. The active components will be integrated in with a small step size (microstep), while the latent components will be integrated with a much larger step size (macrostep).
In summary, this powerful technique (envelope transient over shooting with MRK) can be seen as a multirate scheme (different timestep integration sizes to state variables that present significantly disparate rates of change) coupled with a multirate excitation regime (multiple timescale representations). Consequently, it is able to benefit from the circuits’ heterogeneities, and also from the stimuli timerate disparities, to significantly reduce the computational effort required for simulating the circuits.
3.5. Performance of the Methods
The performance and the efficiency of all the numerical algorithms based on MRK schemes described in this paper were already successfully attested through their application to several illustrative examples of practical relevance. Indeed, electronic circuits with distinct configurations and levels of complexity were especially selected to illustrate the significant gains in computation speed that can be achieved when simulating the circuits with these methods.
Until 2006, multirate RungeKutta methods were only used to obtain the numerical solution of univariate initial value problems (transient electronic circuit simulation in the onedimensional time). For example, in [11, 13], or [15], multirate RungeKutta methods were used to obtain, in an efficient way, the numerical solution of two digital electronic circuits: a digital inverter chain and a pulse generator. Significant reductions in the computational effort were obtained with the methods when compared to standard RK schemes.
Since 2007, multirate RungeKutta methods were also included in numerical algorithms conceived for operating in multivariate frameworks, that is, for solving multidimensional problems described by partial differential systems. For example, in [16, 19] multirate RungeKutta methods were used within a bidimensional envelope following technique to efficiently compute the numerical solution of two RF circuits operating in two distinct time scales: a resistive field effect transistor (FET) mixer and a polar power amplifier used in wireless transmitters. Later, in [20, 21] multirate RungeKutta methods were also used within a threedimensional envelope following technique to simulate a polar power amplifier operating in three separated time scales. The efficiency gains provided by the use of different step sizes (MRK schemes) within the classical univariate timestep integration, but mostly within the bivariate formulation, were also evidenced in [22].
In order to provide a realistic idea of the efficiency (in terms of computational speed) of the methods reviewed in this paper, we decided to include in this section a brief comparison between standard RK and multirate MRK schemes, within the univariate and bivariate formulations. For that, we considered the RF polar transmitter PA described in [16], and depicted in Figure 6, as our illustrative application example. The most relevant components’ values of this circuit are = 25 V, = 2 μH, = 3.2 nF, = 40 nH, = 16 pF, = 0.4 nH, = 16 pF, and R = 50 Ω. The MOSFETs (metaloxide semiconductor field effect transistors) are represented by the following simplified nonlinear device model: with V^{−1}, β=0.25 A/V, , and V. Similarly, the diode currentvoltage characteristic is given by where μ, , and V_{Temp} = 0.026 V.
In this circuit, we have a combination of periodic (RF carrier and digital clock) and aperiodic [AM(t) and PM(t)] forcing functions, of very distinct time scales, with a mixture of heterogeneous state variables with widely disparate rates of variation. For instance, while voltages and currents in the output bandpass filter are all very fast (active state variables), voltages and currents in the AM branch are much slower (latent state variables). The circuit was simulated in MATLAB, and numerical computation times for two different simulation time intervals are presented in Table 1. These results were obtained using (i) a timestep integrator based on a RK scheme (SPICElike simulation), (ii) a timestep integrator based on an MRK scheme, (iii) an envelope transient over shooting (ETS) algorithm using an RK scheme and (iv) an ETS algorithm using an MRK scheme. Standard RK and modern MRK schemes of order 3 based on the BogackiShampine formulas [26] were used to perform the necessary numerical timestep integration.

Table 1 evidences the efficiency gains that can be achieved with the use of MRK schemes, in comparison to classical RK ones, within the univariate and the bivariate formulations. It also shows the advantage of operating with multiple time variables, instead of working in the natural onedimensional time. It must be noted that these gains in computational speed were obtained without compromising the accuracy of the results. Indeed, the maximum discrepancy between the solutions obtained with any of the methods under analysis was of the order 10^{−8} for all the circuit’s state variables.
4. Conclusion
Despite the scientific field of numerical simulation of electronic problems has appeared in the 1970s with the advent of SPICE, this subject is still today a hot topic. Indeed, serious difficulties arise when we have heterogeneous circuits composed of different types of blocks and/or operating in multiple time scales. In this paper, we have briefly reviewed some ways to circumvent these difficulties, operating strictly in the time domain, and which consist in using multirate RungeKutta schemes encapsulated in (i) onedimensional engines and (ii) multidimensional algorithms (by decoupling the components of the signals into different time dimensions). A key aspect that contributes to the success of these numerical methods is the partitioning of the circuits into subcircuits according to the time rates of change of their state variables.
Acknowledgment
The authors would like to acknowledge the financial support provided by the Portuguese Science and Technology Foundation (FCT) under the Project SOPASPTDC/EEATEL/114530/2009.
References
 S. A. Maas, Nonlinear Microwave and RF Circuits, Artech House, Norwood, Mass, USA, 2nd edition, 2003. View at: Zentralblatt MATH
 P. J. Rodrigues, ComputerAided Analysis of Nonlinear Microwave Circuits, Artech House, Norwood, Mass, USA, 1998. View at: Zentralblatt MATH
 E. Hairer, S. P. Nørsett, and G. Wanner, Solving Ordinary Differential Equations. I, vol. 8 of Springer Series in Computational Mathematics, Springer, Berlin, Germany, 1987. View at: Zentralblatt MATH
 E. Hairer and G. Wanner, Solving Ordinary Differential Equations. II, vol. 14 of Springer Series in Computational Mathematics, Springer, Berlin, Germany, 1991. View at: Zentralblatt MATH
 J. D. Lambert, Numerical Methods for Ordinary Differential Systems, John Wiley & Sons Ltd., Chichester, UK, 1991. View at: Zentralblatt MATH
 L. Nagel, SPICE2: A Computer Program to Simulate Semiconductor Circuits, Electronics Research Laboratory, University of California, Berkeley, Calif, USA, 1975. View at: Zentralblatt MATH
 K. Kundert, J. White, and A. SangiovanniVincentelli, SteadyState Methods for Simulating Analog and Microwave Circuits, Kluwer Academic, Norwell, Mass, USA, 1990. View at: Zentralblatt MATH
 J. C. Pedro and N. B. Carvalho, Intermodulation Distortion in Microwave and Wireless Circuits, Artech House, Norwood, Mass, USA, 2003. View at: Zentralblatt MATH
 M. Günther and P. Rentrop, “Multirate ROW methods and latency of electric circuits,” Applied Numerical Mathematics, vol. 13, no. 13, pp. 83–102, 1993. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 M. Günther and P. Rentrop, “Partitioning and multirate strategies in latent electric circuits,” International Series of Numerical Mathematics, vol. 117, pp. 33–60, 1994. View at: Google Scholar  Zentralblatt MATH
 A. Kværnø and P. Rentrop, Low Order Multirate RungeKutta Methods in Electric Circuit Simulation, IWRMM Universität Karlsruhe, 1999.
 A. Kværnø, “Stability of multirate RungeKutta schemes,” International Journal of Differential Equations and Applications, vol. 1A, no. 1, pp. 97–105, 2000. View at: Google Scholar
 M. Günther, A. Kværnø, and P. Rentrop, “Multirate partitioned RungeKutta methods,” BIT, vol. 41, no. 3, pp. 504–514, 2001. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 A. Bartel, M. Günther, and A. Kværnø, “Multirate methods in electrical circuit simulation,” in Progress in Industrial Mathematics at ECMI 2000, pp. 258–265, Springer, 2002. View at: Google Scholar
 J. Oliveira and A. Araújo, “Envelope transient simulation of nonlinear electronic circuits using multirate RungeKutta algorithms,” WSEAS Transactions on Electronics, vol. 3, no. 2, pp. 77–84, 2006. View at: Google Scholar
 J. F. Oliveira and J. C. Pedro, “An efficient timedomain simulation method for multirate RF nonlinear circuits,” IEEE Transactions on Microwave Theory and Techniques, vol. 55, no. 11, pp. 2384–2392, 2007. View at: Publisher Site  Google Scholar
 J. Roychowdhury, “Analyzing circuits with widely separated time scales using numerical PDE methods,” IEEE Transactions on Circuits and Systems, vol. 48, no. 5, pp. 578–594, 2001. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 T. Mei, J. Roychowdhury, T. S. Coffey, S. A. Hutchinson, and D. M. Day, “Robust, stable timedomain methods for solving MPDEs of fast/slow systems,” IEEE Transactions on ComputerAided Design of Integrated Circuits and Systems, vol. 24, no. 2, pp. 226–238, 2005. View at: Publisher Site  Google Scholar
 J. F. Oliveira and J. C. Pedro, “A new timedomain simulation method for highly heterogeneous RF circuits,” Proceedings of the 37th European Microwave Conference (EUMC '07), pp. 1161–1164, 2007. View at: Publisher Site  Google Scholar
 J. F. Oliveira and J. C. Pedro, “An innovative timedomain simulation technique for strongly nonlinear heterogeneous RF circuits operating in diverse timeScales,” in Proceedings of the 38th European Microwave Conference (EuMC '08), pp. 1557–1560, Amsterdam, The Netherlands, October 2008. View at: Publisher Site  Google Scholar
 J. F. Oliveira and J. C. Pedro, “A multipleline double multirate shooting technique for the simulation of heterogeneous RF circuits,” IEEE Transactions on Microwave Theory and Techniques, vol. 57, no. 2, pp. 421–429, 2009. View at: Publisher Site  Google Scholar
 J. F. Oliveira and J. C. Pedro, “Advanced timedomain techniques for strongly nonlinear RF circuit simulation: recent developments and remaining challenges,” in Proceedings of the International Conference on Computer as a Tool (EUROCON '11), Lisbon, Portugal, April 2011. View at: Publisher Site  Google Scholar
 C. E. Christoffersen and J. Alexander, “An adaptive time step control algorithm for nonlinear time domain envelope transient,” in Proceedings of the Canadian Conference on Electrical and Computer Engineering, vol. 2, pp. 883–886, Ontario, Canada, May 2004. View at: Google Scholar
 A. Brambilla and P. Maffezzoni, “Envelopefollowing method to compute steadystate solutions of electrical circuits,” IEEE Transactions on Circuits and Systems, vol. 50, no. 3, pp. 407–417, 2003. View at: Publisher Site  Google Scholar
 E. Ngoya and R. Larcheveque, “Envelop transient analysis: a new method for the transient and steady state analysis of microwave communication circuits and systems,” in IEEE MTTS International Microwave Symposium Digest, vol. 3, pp. 1365–1368, San Francisco, Calif, USA, June 1996. View at: Google Scholar
 P. Bogacki and L. F. Shampine, “A 3(2) pair of RungeKutta formulas,” Applied Mathematics Letters, vol. 2, no. 4, pp. 321–325, 1989. View at: Publisher Site  Google Scholar  Zentralblatt MATH
Copyright
Copyright © 2012 Jorge F. Oliveira and José C. Pedro. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.