Numerical simulation of complex and heterogeneous electronic systems can be a very challenging issue. Circuits composed of a combination of analog, mixed-signal and digital blocks or even radio frequency (RF) blocks, integrated in the same substrate, are very difficult to simulate as a whole at the circuit level. The main reason is because they contain a lot of state variables presenting very distinct properties and evolving in very widely separated time scales. Examples of practical interest are systems-on-a-chip (SoCs), very common in mobile electronics applications, as well as in many other embedded electronic systems. This chapter is intended to briefly review some advanced circuit-level numerical simulation techniques based on circuit-block partitioning schemes, which were especially designed to address the simulation challenges brought by this kind of circuits into the computer-aided-design (CAD) field.
- numerical simulation
- electronic systems
- multirate schemes
- circuit-block partition
Electronic circuit simulation has emerged in the 1970s, triggered by the necessity of engineers having a tool to help in design and analysis of integrated circuits (ICs). Since probing internal nodes of semiconductor chips is extremely difficult, or even prohibitive in almost all cases, manufacturing integrated circuits having no help of a simulation tool would lead to an unbearable set of successive physical prototypes until a final solution was achieved. A simulation package combining device modeling and numerical simulation will help engineers to verify correctness and debug circuits during their design, avoiding physical prototyping and reducing product development expenses. Over the years, the continuous scaling of semiconductor devices and the increasing complexity of electronic architectures have been making CAD tools more and more important for circuits and systems designing. Demands for continuously providing new systems’ functionalities, lower-power consumptions or higher transmission rates (e.g., RF and microwave communication systems) are typical requisites that have led to a complex scenario of highly heterogeneous electronic systems. Conventional algorithms are not capable of simulating such kind of electronic systems in an efficient way. The justification for such ineffectiveness relies on the fact that standard simulation techniques do not perform any distinction between nodes or blocks within the circuits, treating all the variables in the same manner. This causes all the blocks of the circuit (analog, mixed-signal, digital or RF blocks) to be computed with the same numerical scheme, without taking their nature into consideration. To cope with this scenario, some advanced numerical simulation algorithms based on circuit-block partition have been proposed in recent years in the scientific literature. The most important ones are briefly reviewed in this chapter.
2. Review of basic circuit simulation concepts
2.1. Mathematical modeling of electronic systems
Dynamic behavior of electronic systems is modeled as systems of differential algebraic equations (DAEs) involving voltages, currents, charges and fluxes. These systems are usually obtained via nodal analysis, or modified nodal analysis (MNA), which consists of applying the Kirchhoff’s current law (KCL) to each electrical node and writing the branch currents in terms of the circuit node voltages using the corresponding constitutive relations to each circuit element. Such systems have, in general, the following form:
in which is the vector of independent stimuli (voltage or current sources) to the circuit, is the vector of unknowns (voltages and currents waveforms) and
Applying the chain differentiation rule to the reactive term of the DAE system of (1), we are able to obtain
where is habitually denoted as the
which can be rewritten in the classical form
ordinarily utilized in the mathematical literature. When is singular, the DAE system of (3) will not degenerate into a 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 of (5).
2.2. Classic SPICE-like simulation
The most natural way of simulating the dynamic behavior of an electronic circuit is to numerically time-step integrate, in time domain, the DAE system, or the ODE system, modeling its operation. This means that the solution of (1), or (5), has to be computed over a specified time interval from a specific initial condition , leading to the so-called
Computing the solution of (6), or (7), is frequently referred to as
This time-step integration technique was used in the first digital computer programs of circuit analysis initially developed at the Electronics Research Laboratory of the University of California, Berkeley, in the early 1970s, and is still today the most widely used technique for such purpose. It is the core of all SPICE (which means Simulation Program with Integrated Circuit Emphasis) or SPICE-like computer programs, available in many commercial simulators.
2.3. Periodic steady-state simulation
As described earlier, SPICE-like simulation tools are primarily focused on transient analysis. However, in some cases electronics designers are essentially interested in obtaining circuits’ steady-state responses and not their transient regimes. The reason for that is because some specific properties of the circuits are better characterized, or simply only defined, in steady-state (e.g., impedance, voltage gain, current gain, harmonic distortion, signal to noise ratio, etc.).
SPICE tools are not adequate for computing steady-state responses of circuits presenting very different time constants, or high
Periodic steady-state response of an electronic circuit is a regime where its unknowns are a set of generic waveforms (node voltages and branch currents) presenting a common period. Computing the periodic steady-state response, without having to first integrate all the transients, consists of finding the initial condition, , for the DAE, or ODE, system describing the circuit’s operation, such that the values of the unknowns at the end of one period
or, in the ODE form, as
in which is known as the
The most widely used techniques for computing the periodic steady-state solution of electronic circuits are briefly reviewed in the following: the shooting-Newton method  and the harmonic balance method [5, 6].
2.3.1. The shooting-Newton method
The shooting-Newton method  is the time-domain technique most commonly used by the electronic design automation (EDA) community for computing periodic steady-state solutions of electronic circuits. Shooting-Newton is an iterative technique that: (1) starts by computing the solution of the circuit for one period
Shooting is an iterative solver that uses an initial value technique to solve a boundary value problem. With the purpose of providing some technical details on the implementation of the shooting-Newton method, let us consider (8). Since with shooting, we perform numerical time-step integration of the DAE system from until , the main difficulty is that we do not know a priori for which initial condition has to be considered that will lead to the steady-state solution, that is, that will satisfy the periodic boundary condition . Thus, the key aspect of shooting-Newton relies on finding the solution of
Let us now define , where is known as the
Although electronic circuits may operate in strongly nonlinear regimes, their state-transition functions are often moderately nonlinear (or even quite linear). This means that slight perturbations on the initial condition (starting state) produce almost proportional perturbations in the subsequent time states. Taking this aspect into account, it is straightforward to conclude that (11) can be iteratively solved in an efficient way with the Newton’s method, which will lead to the following iterative scheme:
where is the identity matrix. The cost of the solution of the linear system of (12) is dominated by the computational effort required to evaluate the derivative of the state-transition function (usually referred to as the
Although solving (12) and computing the sensitivity matrix (13) involve substantial computational cost, shooting-Newton converges to the steady-state solution much faster than classic time-step integration. This is the reason why it is the time-domain steady-state engine most widely used in the circuit simulation field.
2.3.2. The harmonic balance method
The harmonic balance (HB) method [4, 5, 6] is a frequency domain steady-state simulation technique widely used by the EDA community. In contrast to time domain tools, which represent waveforms as a set of time samples, frequency domain techniques represent periodic signals using coefficients in a sum of complex exponentials (or sines and cosines) harmonically related. The main advantage of HB over time-domain techniques (e.g., shooting-Newton) is that it can represent steady-state solutions (voltage and current waveforms) very accurately using a small number of coefficients. This is especially evident for moderately nonlinear circuits excited by smooth waveforms, in which significant reductions in the computational cost are achieved when HB is used as the simulation tool. However, it must be noted that HB is not suitable for dealing with strongly nonlinear regimes producing waveforms with sharps transitions. In such case, the large number of terms required in the Fourier series expansions will make HB very inefficient.
With the purpose of providing a brief and intuitive explanation of the HB method, let us consider (8). For achieving simplicity in our formulation, let us consider a very small circuit driven by a single periodic source satisfying and whose dynamic behavior is described by a unique unknown, (the generalization to is straightforward). Given that the steady-state regime of the circuit will be periodic with the same period , both the stimulus and the steady-state solution can be represented as the Fourier series
in which is the fundamental frequency. If we substitute (14) into (8), and consider an appropriate harmonic truncation at some order , we will obtain
The HB method converts the differential problem of (8) into the frequency domain, obtaining the algebraic equations system
where , , .
and , are vectors containing the Fourier coefficients of and , respectively. The algebraic equations system of (16) is usually denoted as the
is the composite conversion matrix, known as the
3. Univariate time-domain partitioned simulation engines
3.1. Time-domain latency
As highlighted in the Introduction of this chapter, dynamic regimes of operation of some electronic systems may involve signals (voltages and currents) presenting widely distinct time evolution rates. Typical examples of that are coupled analog-digital systems or combined technologies of baseband analog, digital and RF blocks, in the same circuit. In these examples, very fast signals and slowly varying signals cohabit in the same framework. This property, the one of having in the same problem signals presenting rapid time rates of change, while others evolve in a very slow way (or remain approximately constant within a certain time window) is usually denoted as
In Section 2, we have seen that SPICE-like simulation engines (which are based on time-step integration schemes) are widely used for computing the numerical solution of electronic systems. However, when dealing with circuits presenting time-domain latency, that is, containing node voltages and brunch currents evolving at very different rates, traditional SPICE simulators become inefficient because they expend unnecessary work on the computation of the slowly varying components. This is so because classic initial value solvers (as LMM or RK schemes) integrate all the DAE, or ODE, unknowns with the same step size.
3.2. Partitioned algorithms for time-step integration
To deal with the earlier described time-domain latency in an effective way, some modern partitioned algorithms for univariate time-step integration have been proposed in the recent years in the scientific literature [8, 9, 10]. These powerful techniques, denoted as multirate Runge-Kutta (MRK) schemes, split the ODE system of (5) into coupled fast and slow (latent) subsystems, obtaining
where is the vector containing the fast-varying signals and is the vector holding the slowly varying (latent) ones. The former will be integrated with a small step length (
Consider two Runge-Kutta methods of and stages represented by their Butcher tableaus  and ,
The MRK method conceived for efficiently computing the numerical solution of the partitioned differential system of (19), using a microstep for integrating the fast unknowns and a macrostep for integrating the latent unknowns, is defined as follows [8, 9].
The fast-varying vector is obtained by
The slowly varying vector is given by
From this definition we can attest that numerical coupling between fast and slow differential subsystems is achieved by and . Effective stratagems for computing these intermediate stage values are proposed in [8, 9].
3.3. Circuit-block partitioning strategy
We must be aware that circuit-block partition into fast and slow subsystems can possibly change during the time-step integration process. Hence, it will be very helpful if the simulator is capable of automatically detect the slow and fast variables in the circuit. This automatic classification can be achieved using embedded RK methods  and error estimates usually evaluated for step size control and stiffness detection . Fast-slow partitioning strategies, step size control tools, number of microsteps within a macrostep, stiffness detection stratagems and many other technical details of MRK code implementation are thoroughly addressed in [8, 9, 10] and also in .
Finally, it must be pointed out that significant efficiency gains in computation speed for the simulation of several illustrative examples have been reported in the scientific literature, which demonstrate an effective reduction on the MRK computational cost in comparison to traditional SPICE simulation engines.
4. Multitime partitioned simulation engines
This section is devoted to briefly review some advanced circuit-block partitioning numerical simulation techniques operating in a multivariate time-domain framework. Section 4.1 introduces the multivariate formulation theory, in which the 1-D time is converted into a set of artificial time variables. Section 4.2 addresses some fundamental aspects of the numerical simulation algorithms. Finally, Section 4.3 describes some techniques for automatic circuit-block partition.
4.1. Multivariate formulation
The multivariate formulation is a useful stratagem that plays an important role in the EDA scientific community, especially in the RF and microwave areas. It was initially introduced in 1996  as a sophisticated derivation of quasi-periodic HB, and it has been adopted by other researchers (e.g., [13, 14, 15, 16, 17, 18, 19]), which have demonstrated that it can be an efficient strategy when dealing with electronic circuits operating on widely distinct time scales. The success of multivariate formulation relies on the fact that voltages and currents containing components that evolve themselves at two, or more, widely separated rates of variation can be represented much more efficiently if we define them as functions of two, or more, time variables (artificial time scales). With this stratagem all signals (stimuli and responses) will be represented as multivariate functions, which will imply that dynamic behavior of the circuits will no longer be modeled by DAE, or ODE, systems formulated in the 1-D time. It will be described by partial differential systems.
In order to provide a simple and illustrative mathematical description of the multivariate formulation let us consider the classical example of a generic nonlinear RF circuit driven by an envelope modulated signal of the form
where and are, respectively, the amplitude and phase slowly varying baseband signals (envelope), modulating the fast-varying carrier. Computing the numerical solution of such circuit using conventional 1-D time-step integrators (RK methods or LMM methods) is very expensive. This is so because the response of the circuit has to be evaluated during a long time window defined by the slowly varying envelope, wherein the time-step size is severely restricted by the high frequency carrier. Taking into consideration the disparity between the baseband and the carrier time scales, and assuming that they are also uncorrelated, which is typically true, we are able to rewrite (26) as a bivariate function
where is the slow baseband artificial time scale and is the fast carrier artificial time scale. It must be noted that is a periodic function with respect to but not to , that is,
This means that a generic 1-D time interval will be mapped into a 2-D rectangular domain. It is easy to attest that, in general, the number of points required to represent in is much less than the number of points needed to represent in . This is especially evident when the and time scales are widely separated [13, 14].
Let us now consider the DAE system of (1) modeling the dynamic behavior of a generic RF circuit excited by the envelope modulated signal of (26). Taking the abovementioned into account, we will adopt the following procedure for the vector-valued functions and : is replaced by in the slowly varying entities (envelope time scale) and is replaced by in the fast-varying entities (RF carrier time scale). The application of this stratagem will turn the DAE system of (1) into the following partial differential system
usually denoted as
where represents the remainder of division of by .
The generalization of the bivariate strategy to a multidimensional problem with more than two time scales is straightforward. In fact, if the signals in the circuit present separate rates of change, then time scales will be used. In that case (29) assumes the generic form
and the univariate solution, , may be recovered from its multivariate form, , by setting .
4.2. Partitioned algorithms for envelope following computation
The main advantage of the earlier described MPDAE approach is that it can result in significant improvements in simulation speed when compared to DAE-based alternatives [13, 14, 15]. However, by itself this approach does not perform any distinction between nodes or blocks in the circuit under analysis. In fact, in the first multivariate circuit simulation schemes initially introduced in , and then in , the same numerical algorithm was used to compute all the unknowns of the circuit. Only a few years later other advanced multivariate algorithms were proposed (e.g., [16, 17, 18, 19]) in way to take into account possible circuit’s heterogeneities. These algorithms are based on circuit-block partitioning stratagems defined within the multivariate frameworks. The most important ones regarding pure time-domain operations are briefly reviewed in the following.
As stated earlier, envelope modulated regimes are typical cases of practical interest. Computing responses to excitations of the form of (26) is a technique generally referred to as
defined on the rectangular domain . is a vector-valued initial-condition satisfying , and is the periodic boundary condition due to the periodicity of the problem in the fast time scale. In order to solve this initial-boundary value problem let us begin to consider the following semi-discretization of in the slow time dimension defined by
in which is the total number of steps in . Now, using a backward differentiation formula (BDF)  to approximate the derivatives of the MPDAE (for simplicity let us consider here the Gear-2 rule ), we obtain for each slow time instant , from to the periodic boundary value problem defined by
where . This means that, once is known, the solution on the next slow time instant, , is obtained by solving (34). Hence, a set of boundary value problems have to be solved for obtaining the whole bivariate solution in the entire domain . A very powerful technique has been proposed in the literature for solving each of the periodic boundary problems defined by (34) in an efficient way [16, 17]. This technique uses shooting-Newton based on MRK schemes. It splits the circuits into two distinct subsets according to the time evolution rates of their voltages and currents and performs time-step integration with different step lengths in each of the consecutive shooting iterations needed to solve (34). For that (34) is firstly divided into the following coupled fast-slow subsystems
is the vector containing the fast-varying circuit variables at the slow time instant , and is the vector holding the slowly varying ones at the same slow time instant. The former will be integrated in with a small step size
4.3. Circuit-block partitioning strategy
Since the partition of the electronic system in fast and slow subsystems may dynamically vary with time, it will be of great usefulness if the simulation algorithm is capable of automatically detecting the fast-varying and the slowly varying signals. The approach adopted in [16, 17] for that purpose is briefly described in the following.
As mentioned earlier, each of the periodic boundary value problems defined by (34) is solved using shooting-Newton based on multirate time-step integrators (MRK schemes). Now, taking into account that shooting is an iterative technique, the key idea of this partitioning strategy is to use a uni-rate scheme on the first shooting iterations needed for each time line of the rectangular domain. For that, it starts by considering in the adopted MRK scheme, so that it degenerates into a standard uni-rate RK scheme. This means that all the unknowns (voltages or currents) will be integrated in with the same microstep . After that, the differential system of (34) is partitioned into (35) according to the variations in the derivatives of the unknowns. Each unknown that practically evidences no variations in its time derivatives for the entire time line , that is, each unknown that satisfies the condition
where and is a small specified tolerance, will be treated as slow on the next shooting iterations. The remaining unknowns will be treated as fast. In a cheap scheme within a uniform grid each slope can be simply given by
All the subsequent shooting iterations on the time line will be conducted in a multirate way using different step sizes. In a nonuniform grid, step sizes and may be chosen using any step-size control tool. In a uniform grid, they can be predefined or successively refined to achieve a desired accuracy.
The robustness of this partitioning strategy can be improved if more than one shooting iteration with the uni-rate scheme is considered for each time line . However, such tactic will obviously conduct to some extra computational cost, turning the overall algorithm to be a little less efficient.
Finally, to conclude this section, it must be highlighted that although efficiency gains are essentially dependent on the ratio between the number of slow and fast signals, significant reductions on the computational effort were reported in the scientific literature (e.g., [16, 17]) for the simulation of several illustrative examples of practical interest using the envelope following bivariate partitioned techniques.
5. Hybrid time-frequency partitioned techniques
This section is intended to provide a brief description of powerful circuit-block partitioning numerical simulation techniques operating in multivariate hybrid time-frequency frameworks. Section 5.1 introduces the fundaments of multivariate hybrid envelope following. Section 5.2 addresses some basic details of the simulation algorithms and finally, Section 5.3 describes the approach for automatic circuit-block partition.
5.1. Multivariate hybrid envelope following
Let us once again consider the initial-boundary value problem of (32) characterizing the bivariate nature of an electronic circuit operating at two widely separated time scales. Let us also consider the semi-discretization of the rectangular domain, which has led to the set of periodic boundary value problems defined by (34). When waveforms are not excessively demanding on the quantity of harmonic components needed for an accurate frequency domain representation, we are able to use the HB method for efficiently computing the solution of (34). In such case we will obtain for each time line the following HB system
in which is a vector containing the Fourier coefficients of the stimuli (independent sources) and a vector containing the Fourier coefficients of the unknowns (node voltages and brunch currents’ waveforms), at . is the maximum harmonic order and is the number of unknowns. In order to obtain the solution in the entire rectangular domain a total of HB systems have to be solved.
With this hybrid time-frequency approach we have an envelope following technique that handles the slow variations of the unknowns in the time domain and their fast variations in the frequency domain. This technique, usually denoted as multivariate
5.2. Partitioned algorithms for hybrid envelope following computation
The hybrid time-frequency envelope following technique presented earlier does not perform any distinction between nodes or blocks in the circuit. Although it has been widely used, especially by the RF and microwave community, only a few years ago other versions regarding circuit-block partition were proposed [18, 19] in way to take into account possible circuit’s heterogeneities. The most important aspects of these partitioning techniques are briefly reviewed in the following.
Let us rewrite (39) as
The Newton iterative solver is usually utilized to solve (40), leading to
which implies that for computing from the previous estimate we have to solve a linear system composed of equations. The system of (41) involves the so-called Jacobian matrix of , which has a block structure, consisting of an matrix of square submatrices (blocks), each one with dimensions . The general block of row and column can be expressed as
where and denote the Toeplitz matrices  of the vectors containing the Fourier coefficients of and , respectively.
A very powerful technique has been proposed in the literature [18, 19] for solving (41) in an efficient way. This technique takes into account that, in some cases (e.g., RF heterogeneous systems), there are parts of the circuits in which there are no fluctuations dictated by the fast carrier. As a consequence, bidimensional forms of voltages and currents in those parts have no dependence on . This means that for each time line each signal has a constant value that can be represented as a Fourier series with a unique coefficient. Thus, while quickly varying signals are represented in the frequency domain by a set of Fourier coefficients, signals in latent (slowly varying) blocks are represented by a single coefficient. Having this in mind it is straightforward to conclude that the size of vector will be significantly reduced, as well as the size (number of equations) of the HB system. Significant Jacobian matrix size reductions will also be achieved, since some of the submatrices will no longer be matrices, but, instead, simple scalar elements. For instance, let us suppose that the component of (40) is a slowly varying signal exclusively dependent on other slowly varying signals. Let us also suppose that the component is also a slowly varying signal. Then, the Jacobian matrix block of row and column will be a block given by
It may be noted that this scalar block can be seen as a particular case of the general block of (42), if is assumed as the maximum harmonic order. Actually, given that is a constant function evidencing no fluctuations in the fast time scale, the last term of (42) will vanish and there will be no more necessity of converting the right-hand side terms of (43) into the frequency domain.
Since only fast-varying signals are converted into the frequency domain, this partitioned technique can be seen as a combination of multivariate envelope transient harmonic balance (by treating the fast-varying signals in a hybrid time-frequency framework) with a pure time marching simulation engine (by treating some of the signals in a pure time domain scheme). Thus, beyond the notorious significant vector and matrix size reductions above mentioned, there is another important advantage brought by this partitioned technique. For example, in complex heterogeneous RF systems strongly nonlinear regimes of operation are in general associated to digital or baseband blocks, whereas moderately nonlinear regimes are typical of RF blocks. With this partitioned technique signals in digital and baseband blocks are appropriately computed in the time domain, while signals in RF blocks are treated in the hybrid time-frequency framework.
5.3. Circuit-block partitioning strategy
Similar to what we have mentioned for the methods discussed in the previous sections, it will be of great utility if the simulator is able to automatically distinguish the fast-varying variables from the slowly varying ones. We now briefly review the approach addressed in  for that purpose, which splits the circuit into distinct blocks according to the time rates of change of their voltages and currents.
Let us consider the HB system of (39). As stated earlier, this system has to be solved with the iterative scheme of (41) for each artificial time line . The automatic partitioning strategy proposed in  consists in considering all the circuit’s variables as fast on the first iteration of (41), that is, it consists in initially representing all the unknowns in the frequency domain as a set of Fourier coefficients. This way, the algorithm starts by computing a single iteration in (41) to evaluate
where each , , is a vector defined as
Each of the is then inspected. If its Fourier coefficients of order are practically null (their absolute values stay under a very small prescribed tolerance) it will be classified as slow. Otherwise it will be classified as fast. After this classification, which will temporarily split the system into fast and slow subsystems, the simulator considers that signals in slow subsystems can be represented by a single Fourier coefficient for the remaining iterations needed to compute the solution of (41).
Similar to what we have discussed for the method described in Section 4, the robustness of this partitioning strategy may be improved if more than one iteration in (41) is computed before the simulator decides which signals are slow and which signals are fast. The main drawback of such approach is the loss of some efficiency due to the extra computational effort required.
As a final point of this section, we would like to point out that significant gains in computation speed have been reported in the scientific literature (e.g., [18, 19]) for the simulation of several illustrative examples of practical interest using these hybrid time-frequency envelope following partitioned algorithms.
In this chapter, we have briefly reviewed some powerful numerical simulation techniques based on partitioned stratagems. Such techniques were especially designed to cope with the simulation challenges brought by emerging electronic technologies to the EDA community, as is the case of complex heterogeneous electronic systems composed of a combination of different kinds of circuit blocks (analog, mixed-signal and digital blocks, or even radio frequency blocks) containing node voltages and brunch currents of very distinct formats and running on widely separated time scales. With these partitioned techniques signals within different blocks of the circuits are computed with distinct algorithms and/or step sizes. Considerable reductions on the computational cost have been registered in several experiments published in the scientific literature (in comparison to previously recognized techniques) without compromising the accuracy of the results.
This work is funded by National Funds through FCT—Fundação para a Ciência e Tecnologia, under the project UID/EEA/50008/2013.
Lambert J. Numerical Methods for Ordinary Differential Systems: The Initial Value Problem. 1st ed. West Sussex: Wiley; 1991
Hairer E, Nørsett S, Wanner G. Solving Ordinary Differential Equations I: Nonstiff Problems. 1st ed. Berlin: Springer-Verlag; 1987. DOI: 10.1007/978-3-662-12607-3
Hairer E, Wanner G. Solving Ordinary Differential Equations II: Stiff and Differential Algebraic Problems. 2nd ed. Berlin: Springer-Verlag; 1996. DOI: 10.1007/978-3-642-05221-7
Kundert K, White J, Sangiovanni-Vincentelli A. Steady-State Methods for Simulating Analog and Microwave Circuits. Norwell, MA: Spinger; 1990
Rodrigues P. Computer-Aided Analysis of Nonlinear Microwave Circuits. Norwood, MA: Artech House; 1998
Maas S. Nonlinear Microwave and RF Circuits. 2nd ed. Norwood, MA: Artech House; 2003
Pedro J, Carvalho N. Intermodulation Distortion in Microwave and Wireless Circuits. 1st ed. Norwood MA: Artech House; 2003
Kværnø A. Stability of multirate Runge-Kutta schemes. International Journal of Differential Equations and Applications. 2000; 1A:97-105
Günther M, Kværnø A, Rentrop P. Multirate partitioned Runge-Kutta methods. BIT. 2001; 41(3):504-514. DOI: 10.1023/A:1021967112503
Bartel A, Günther M, Kværnø A. Multirate methods in electrical circuit simulation. Progress in Industrial Mathematics, Springer. 2002; 1:258-265. DOI: 10.1007/978-3-662-04784-2_35
Oliveira J, Araújo A. Envelope transient simulation of nonlinear electronic circuits using multirate Runge-Kutta algorithms. WSEAS Transactions on Electronics. 2006; 3(2):77-84
Brachtendorf H, Welsch G, Laur R, Bunse-Gerstner A. Numerical steady-state analysis of electronic circuits driven by multi-tone signals. Electrical Engineering (Springer-Velag). 1996; 79(2):103-112. DOI: 10.1007/BF01232919
Roychowdhury J. Analyzing circuits with widely separated time scales using numerical PDE methods. IEEE Transactions on Circuits and Systems. 2001; 48(5):578-594. DOI: 10.1109/81.922460
Mei T, Roychowdhury J, Coffey T, Hutchinson S, Day D. Robust stable time-domain methods for solving MPDEs of fast/slow systems. IEEE Trans. on Computer-Aided Design of Integrated Circuits and Systems. 2005; 24(2):226-239. DOI: 10.1109/TCAD.2004.841073
Zhu L, Christoffersen C. Adaptive harmonic balance analysis of oscillators using multiple time scales. In: 3rd International IEEE Northeast Workshop on Circuits and Systems; Québec City. 2005. pp. 187-190
Oliveira J, Pedro J. An efficient time-domain simulation method for multirate RF nonlinear circuits. IEEE Transactions on Microwave Theory and Techniques. 2007; 55(11):2384-2392. DOI: 10.1109/TMTT.2007.908679
Oliveira J, Pedro J. A multiple-line double multirate shooting technique for the simulation of heterogeneous RF circuits. IEEE Transactions on Microwave Theory and Techniques. 2009; 57(2):421-429. DOI: 10.1109/TMTT.2008.2011228
Oliveira J, Pedro J. A new mixed time-frequency simulation method for nonlinear heterogeneous multirate RF circuits. In: Proceeding of the IEEE MTT-S International Microwave Symposium (MTT ‘10); May 2010; Anaheim, CA. pp. 548-551
Oliveira J, Pedro J. Efficient RF circuit simulation using an innovative mixed time-frequency method. IEEE Transactions on Microwave Theory and Techniques. 2011; 59(4):827-836. DOI: 10.1109/TMTT.2010.2095035
Rizzoli V, Neri A, Mastri F. A modulation-oriented piecewise harmonic-balance technique suitable for transient analysis and digitally modulated analysis. In: Proc. 26th European Microwave Conference; Jun. 1996; Prague. 1996. pp. 546-550
Sharrit D. Method for simulating a circuit. U.S. Patent 5588142, December 24, 1996