1. Introduction

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

© 2016 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and eproduction in any medium, provided the original work is properly cited. © 2018 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 circuitblock 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:

$$p(y(t)) + \frac{d\boldsymbol{q}(y(t))}{dt} = \mathbf{x}(t),\tag{1}$$

M y½ � ð Þt

dyð Þt

which can be rewritten in the classical form

2.2. Classic SPICE-like simulation

p yð Þþ ð Þt

dyð Þt

(accuracy) and numerical stability.

dq yð Þ ð Þt

problems

or

to convert (3) into the following ordinary differential equations' (ODE) system,

dt <sup>¼</sup> M y½ � ð Þ<sup>t</sup> �<sup>1</sup>

dyð Þt

algebraic equations combined with a set of differential equations of the form of (5).

dyð Þt

where M y½ � ð Þt is habitually denoted as the mass matrix. If M y½ � ð Þt is invertible then it is possible

ordinarily utilized in the mathematical literature. When M y½ � ð Þt 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

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 t0; tFinal ½ � from a specific initial condition yð Þ¼ t<sup>0</sup> y0, leading to the so-called initial value

Computing the solution of (6), or (7), is frequently referred to as transient analysis and can be done by using initial value solvers, such as linear multistep methods (LMM) [1] or one-step methods (the popular Runge-Kutta (RK) schemes) [2, 3]. Either LMM or RK methods offer a large variety of explicit and implicit schemes, with very distinct properties in terms of order

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.

dt <sup>¼</sup> <sup>x</sup>ð Þ� <sup>t</sup> p yð Þ ð Þ<sup>t</sup> , (3)

http://dx.doi.org/10.5772/intechopen.75490

25

Numerical Simulation of Electronic Systems Based on Circuit-Block Partitioning Strategies

ð Þ xð Þ� t p y½ � ð Þt , (4)

dt <sup>¼</sup> <sup>f</sup>ð Þ <sup>t</sup>; <sup>y</sup>ð Þ<sup>t</sup> , (5)

dt <sup>¼</sup> <sup>x</sup>ð Þ<sup>t</sup> , <sup>y</sup>ð Þ¼ <sup>t</sup><sup>0</sup> <sup>y</sup>0, t<sup>0</sup> <sup>≤</sup> <sup>t</sup> <sup>≤</sup> tFinal, <sup>y</sup>ð Þ<sup>t</sup> <sup>∈</sup> <sup>R</sup>n, (6)

dt <sup>¼</sup> <sup>f</sup>ð Þ <sup>t</sup>; <sup>y</sup>ð Þ<sup>t</sup> , <sup>y</sup>ð Þ¼ <sup>t</sup><sup>0</sup> <sup>y</sup>0, t<sup>0</sup> <sup>≤</sup> <sup>t</sup> <sup>≤</sup> tFinal, <sup>y</sup>ð Þ<sup>t</sup> <sup>∈</sup> <sup>R</sup><sup>n</sup>: (7)

in which <sup>x</sup>ð Þ<sup>t</sup> <sup>∈</sup> <sup>R</sup><sup>n</sup> is the vector of independent stimuli (voltage or current sources) to the circuit, <sup>y</sup>ð Þ<sup>t</sup> <sup>∈</sup> <sup>R</sup><sup>n</sup> is the vector of unknowns (voltages and currents waveforms) and <sup>n</sup> is the total number of unknowns. <sup>p</sup> : <sup>R</sup><sup>n</sup> ! <sup>R</sup><sup>n</sup> describes the memoryless elements in the circuit (linear and nonlinear elements, as resistors, nonlinear voltage-controlled current sources, etc.), while <sup>q</sup> : <sup>R</sup><sup>n</sup> ! <sup>R</sup><sup>n</sup> models all linear and nonlinear reactive circuit elements, as capacitors or inductors, represented as voltage-dependent electric charges or current-dependent magnetic fluxes, respectively.

Applying the chain differentiation rule to the reactive term of the DAE system of (1), we are able to obtain

$$\frac{dq(y)}{dy}\frac{dy(t)}{dt} = \mathfrak{x}(t) - \mathfrak{p}(\mathfrak{y}(t)),\tag{2}$$

or,

Numerical Simulation of Electronic Systems Based on Circuit-Block Partitioning Strategies http://dx.doi.org/10.5772/intechopen.75490 25

$$\mathbf{M}[y(t)]\frac{dy(t)}{dt} = \mathbf{x}(t) - \mathbf{p}(y(t)),\tag{3}$$

where M y½ � ð Þt is habitually denoted as the mass matrix. If M y½ � ð Þt is invertible then it is possible to convert (3) into the following ordinary differential equations' (ODE) system,

$$\frac{dy(t)}{dt} = \mathbf{M}[y(t)]^{-1}(\mathbf{x}(t) - \mathbf{p}[y(t)]).\tag{4}$$

which can be rewritten in the classical form

$$\frac{dy(t)}{dt} = f(t, y(t)),\tag{5}$$

ordinarily utilized in the mathematical literature. When M y½ � ð Þt 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 t0; tFinal ½ � from a specific initial condition yð Þ¼ t<sup>0</sup> y0, leading to the so-called initial value problems

$$p(\mathbf{y}(t)) + \frac{d\mathbf{q}(\mathbf{y}(t))}{dt} = \mathbf{x}(t), \quad \mathbf{y}(t\_0) = \mathbf{y}\_{0^\prime} \quad t\_0 \le t \le t\_{\text{Final}}, \quad \mathbf{y}(t) \in \mathbb{R}^n,\tag{6}$$

or

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 circuitblock partition have been proposed in recent years in the scientific literature. The most impor-

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

dq yð Þ ð Þt

in which <sup>x</sup>ð Þ<sup>t</sup> <sup>∈</sup> <sup>R</sup><sup>n</sup> is the vector of independent stimuli (voltage or current sources) to the circuit, <sup>y</sup>ð Þ<sup>t</sup> <sup>∈</sup> <sup>R</sup><sup>n</sup> is the vector of unknowns (voltages and currents waveforms) and <sup>n</sup> is the total number of unknowns. <sup>p</sup> : <sup>R</sup><sup>n</sup> ! <sup>R</sup><sup>n</sup> describes the memoryless elements in the circuit (linear and nonlinear elements, as resistors, nonlinear voltage-controlled current sources, etc.), while <sup>q</sup> : <sup>R</sup><sup>n</sup> ! <sup>R</sup><sup>n</sup> models all linear and nonlinear reactive circuit elements, as capacitors or inductors, represented as voltage-dependent electric charges or current-dependent magnetic

Applying the chain differentiation rule to the reactive term of the DAE system of (1), we are

dt <sup>¼</sup> <sup>x</sup>ð Þ<sup>t</sup> , (1)

dt <sup>¼</sup> <sup>x</sup>ð Þ� <sup>t</sup> p yð Þ ð Þ<sup>t</sup> , (2)

tant ones are briefly reviewed in this chapter.

24 Numerical Simulations in Engineering and Science

2. Review of basic circuit simulation concepts

element. Such systems have, in general, the following form:

p yð Þþ ð Þt

dq yð Þ dy

dyð Þt

2.1. Mathematical modeling of electronic systems

fluxes, respectively.

able to obtain

or,

$$\frac{dy(t)}{dt} = f(t, y(t)), \quad y(t\_0) = y\_{0^\prime} \quad t\_0 \le t \le t\_{\text{Final}}, \quad y(t) \in \mathbb{R}^n. \tag{7}$$

Computing the solution of (6), or (7), is frequently referred to as transient analysis and can be done by using initial value solvers, such as linear multistep methods (LMM) [1] or one-step methods (the popular Runge-Kutta (RK) schemes) [2, 3]. Either LMM or RK methods offer a large variety of explicit and implicit schemes, with very distinct properties in terms of order (accuracy) and numerical stability.

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

guessed initial condition (which is in general determined from a previous DC analysis); (2) 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 updated and (3) the circuit is then re-simulated for one period with the adjusted initial condition, and this process is repeated until the solution

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 timestep integration of the DAE system from t ¼ t<sup>0</sup> until t ¼ t<sup>0</sup> þ T, the main difficulty is that we do not know a priori for which initial condition yð Þ t<sup>0</sup> has to be considered that will lead to the steady-state solution, that is, that will satisfy the periodic boundary condition yð Þ¼ t<sup>0</sup>

Let us now define ϕð Þ¼ yð Þ t<sup>0</sup> ; T yð Þ t<sup>0</sup> þ T , where ϕ is known as the state-transition function [4, 7],

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

� I

where I is the n � n 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 statetransition function (usually referred to as the sensitivity matrix). This matrix is computed taking into account the chain differentiation rule, that is, taking into consideration that ϕð Þ yð Þ t<sup>0</sup> ; T is, in fact, the numerical vector yK, with K being the total number of time steps in the interval ½ � <sup>t</sup>0; <sup>t</sup><sup>0</sup> <sup>þ</sup> <sup>T</sup> , which depends on the previous step value <sup>y</sup><sup>K</sup>�<sup>1</sup>, which, itself, depends on <sup>y</sup><sup>K</sup>�<sup>2</sup>, and

> <sup>¼</sup> <sup>∂</sup>y<sup>K</sup> <sup>∂</sup>y<sup>K</sup>�<sup>1</sup> � <sup>∂</sup>y<sup>K</sup>�<sup>1</sup> <sup>∂</sup>y<sup>K</sup>�<sup>2</sup>

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

� � � �

yð Þ¼ t<sup>0</sup> y½ �<sup>r</sup> ð Þ t<sup>0</sup>

� ⋯ � ∂y<sup>1</sup> ∂y<sup>0</sup>

yð Þ¼ t<sup>0</sup> yð Þ t<sup>0</sup> þ T ⇔ yð Þ� t<sup>0</sup> yð Þ¼ t<sup>0</sup> þ T 0: (10)

Numerical Simulation of Electronic Systems Based on Circuit-Block Partitioning Strategies

http://dx.doi.org/10.5772/intechopen.75490

27

ϕð Þ� yð Þ t<sup>0</sup> ; T yð Þ¼ t<sup>0</sup> 0: (11)

<sup>y</sup>½ � <sup>r</sup>þ<sup>1</sup> ð Þ� <sup>t</sup><sup>0</sup> <sup>y</sup>½ �<sup>r</sup> ð Þ <sup>t</sup><sup>0</sup>

h i <sup>¼</sup> <sup>0</sup>, (12)

: (13)

yð Þ t<sup>0</sup> þ T . Thus, the key aspect of shooting-Newton relies on finding the solution of

<sup>∂</sup><sup>ϕ</sup> <sup>y</sup>½ �<sup>r</sup> ð Þ <sup>t</sup><sup>0</sup> ; <sup>T</sup> � � ∂yð Þ t<sup>0</sup>

∂ϕð Þ yð Þ t<sup>0</sup> ; T ∂yð Þ t<sup>0</sup>

" #�

at the end of the period matches the initial condition.

and rewrite (10) as

<sup>ϕ</sup> <sup>y</sup>½ �<sup>r</sup> ð Þ <sup>t</sup><sup>0</sup> ; <sup>T</sup>

lead to the following iterative scheme:

� � � <sup>y</sup>½ �<sup>r</sup> ð Þþ <sup>t</sup><sup>0</sup>

so forth. The sensitivity matrix is then given by

widely used in the circuit simulation field.

SPICE tools are not adequate for computing steady-state responses of circuits presenting very different time constants, or high Q resonances, as is typically the case of RF and microwave circuits. This is so because they have to pass through the lengthy process of integrating all transients, and expecting them to vanish. Indeed, in such cases time-step integration can be extremely inefficient, since the number of discretization time steps used by the numerical integration scheme will be dramatically enormous. This is 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 takes to achieve steady-state, while the length of the time steps is constrained by the highest frequency component.

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, yð Þ t<sup>0</sup> , for the DAE, or ODE, system describing the circuit's operation, such that the values of the unknowns at the end of one period T match their values at the beginning of that period, that is to say, yð Þ¼ t<sup>0</sup> yð Þ t<sup>0</sup> þ T . These problems (evaluating the solution to a differential system that satisfies constraints at two or more distinct time instants) are referred to as boundary value problems. In this case, we have a periodic boundary value problem that can be formulated as:

$$p(\mathbf{y}(t)) + \frac{d\mathbf{q}(\mathbf{y}(t))}{dt} = \mathbf{x}(t), \quad \mathbf{y}(t\_0) = \mathbf{y}(t\_0 + T), \quad t\_0 \le t \le t\_0 + T, \quad \mathbf{y}(t) \in \mathbb{R}^n. \tag{8}$$

or, in the ODE form, as

$$\frac{dy(t)}{dt} = f(t, y(t)), \quad y(t\_0) = y(t\_0 + T), \quad t\_0 \le t \le t\_0 + T, \quad y(t) \in \mathbb{R}^n,\tag{9}$$

in which yð Þ¼ t<sup>0</sup> yð Þ t<sup>0</sup> þ T is known as the periodic boundary condition.

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 [4] and the harmonic balance method [5, 6].

### 2.3.1. The shooting-Newton method

The shooting-Newton method [4] 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 T using a LMM or RK integrator, considering some guessed initial condition (which is in general determined from a previous DC analysis); (2) 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 updated and (3) the circuit is then re-simulated for one period with the adjusted initial condition, and this process is repeated until the solution at the end of the period matches the initial condition.

2.3. Periodic steady-state simulation

26 Numerical Simulations in Engineering and Science

constrained by the highest frequency component.

value problem that can be formulated as:

dq yð Þ ð Þt

in which yð Þ¼ t<sup>0</sup> yð Þ t<sup>0</sup> þ T is known as the periodic boundary condition.

p yð Þþ ð Þt

dyð Þt

monic balance method [5, 6].

2.3.1. The shooting-Newton method

or, in the ODE form, as

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 Q resonances, as is typically the case of RF and microwave circuits. This is so because they have to pass through the lengthy process of integrating all transients, and expecting them to vanish. Indeed, in such cases time-step integration can be extremely inefficient, since the number of discretization time steps used by the numerical integration scheme will be dramatically enormous. This is 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 takes to achieve steady-state, while the length of the time steps is

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, yð Þ t<sup>0</sup> , for the DAE, or ODE, system describing the circuit's operation, such that the values of the unknowns at the end of one period T match their values at the beginning of that period, that is to say, yð Þ¼ t<sup>0</sup> yð Þ t<sup>0</sup> þ T . These problems (evaluating the solution to a differential system that satisfies constraints at two or more distinct time instants) are referred to as boundary value problems. In this case, we have a periodic boundary

dt <sup>¼</sup> <sup>x</sup>ð Þ<sup>t</sup> , <sup>y</sup>ð Þ¼ <sup>t</sup><sup>0</sup> <sup>y</sup>ð Þ <sup>t</sup><sup>0</sup> <sup>þ</sup> <sup>T</sup> , t<sup>0</sup> <sup>≤</sup> <sup>t</sup> <sup>≤</sup> <sup>t</sup><sup>0</sup> <sup>þ</sup> T, <sup>y</sup>ð Þ<sup>t</sup> <sup>∈</sup> <sup>R</sup>n, (8)

dt <sup>¼</sup> <sup>f</sup>ð Þ <sup>t</sup>; <sup>y</sup>ð Þ<sup>t</sup> , <sup>y</sup>ð Þ¼ <sup>t</sup><sup>0</sup> <sup>y</sup>ð Þ <sup>t</sup><sup>0</sup> <sup>þ</sup> <sup>T</sup> , t<sup>0</sup> <sup>≤</sup> <sup>t</sup> <sup>≤</sup> <sup>t</sup><sup>0</sup> <sup>þ</sup> T, <sup>y</sup>ð Þ<sup>t</sup> <sup>∈</sup> <sup>R</sup>n, (9)

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 [4] and the har-

The shooting-Newton method [4] 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 T using a LMM or RK integrator, considering some 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 timestep integration of the DAE system from t ¼ t<sup>0</sup> until t ¼ t<sup>0</sup> þ T, the main difficulty is that we do not know a priori for which initial condition yð Þ t<sup>0</sup> has to be considered that will lead to the steady-state solution, that is, that will satisfy the periodic boundary condition yð Þ¼ t<sup>0</sup> yð Þ t<sup>0</sup> þ T . Thus, the key aspect of shooting-Newton relies on finding the solution of

$$\mathbf{y}(t\_0) = \mathbf{y}(t\_0 + T) \Leftrightarrow \mathbf{y}(t\_0) - \mathbf{y}(t\_0 + T) = \mathbf{0}.\tag{10}$$

Let us now define ϕð Þ¼ yð Þ t<sup>0</sup> ; T yð Þ t<sup>0</sup> þ T , where ϕ is known as the state-transition function [4, 7], and rewrite (10) as

$$
\phi(\mathbf{y}(t\_0), T) - \mathbf{y}(t\_0) = \mathbf{0}.\tag{11}
$$

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:

$$\left[\phi\left(\mathbf{y}^{[r]}(t\_0), T\right) - \mathbf{y}^{[r]}(t\_0) + \left[\frac{\partial\phi\left(\mathbf{y}^{[r]}(t\_0), T\right)}{\partial\mathbf{y}(t\_0)} - I\right]\right]\_{\mathbf{y}(t\_0) = \mathbf{y}^{[r]}(t\_0)} \left[\mathbf{y}^{[r+1]}(t\_0) - \mathbf{y}^{[r]}(t\_0)\right] = 0,\tag{12}$$

where I is the n � n 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 statetransition function (usually referred to as the sensitivity matrix). This matrix is computed taking into account the chain differentiation rule, that is, taking into consideration that ϕð Þ yð Þ t<sup>0</sup> ; T is, in fact, the numerical vector yK, with K being the total number of time steps in the interval ½ � <sup>t</sup>0; <sup>t</sup><sup>0</sup> <sup>þ</sup> <sup>T</sup> , which depends on the previous step value <sup>y</sup><sup>K</sup>�<sup>1</sup>, which, itself, depends on <sup>y</sup><sup>K</sup>�<sup>2</sup>, and so forth. The sensitivity matrix is then given by

$$\frac{\partial \phi(\boldsymbol{y}(t\_0), T)}{\partial \boldsymbol{y}(t\_0)} = \frac{\partial \boldsymbol{y}\_K}{\partial \boldsymbol{y}\_{K-1}} \cdot \frac{\partial \boldsymbol{y}\_{K-1}}{\partial \boldsymbol{y}\_{K-2}} \cdot \dots \cdot \frac{\partial \boldsymbol{y}\_1}{\partial \boldsymbol{y}\_0}.\tag{13}$$

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–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 x tð Þ ∈ R satisfying x tðÞ¼ x tð Þ þ T and whose dynamic behavior is described by a unique unknown, y tð Þ∈ R (the generalization to x tð Þ<sup>∈</sup> <sup>R</sup>n, ytð Þ <sup>∈</sup> <sup>R</sup><sup>n</sup> is straightforward). Given that the steady-state regime of the circuit will be periodic with the same period T, both the stimulus and the steady-state solution can be represented as the Fourier series

$$\mathbf{x}(t) = \sum\_{k=-\circ\circ}^{+\circ\circ} \mathbf{X}\_k \mathbf{e}^{j\mathbf{k}\omega\_0 t}, \qquad \mathbf{y}(t) = \sum\_{k=-\circ\circ}^{+\circ\circ} \mathbf{Y}\_k \mathbf{e}^{j\mathbf{k}\omega\_0 t} \tag{14}$$

in which

3.1. Time-domain latency

dF Yð Þ

corresponding to gyt ð Þ¼ ð Þ dp y t ð Þ ð Þ =dy and cyt ð Þ¼ ð Þ dq y t ð Þ ð Þ =dy.

3. Univariate time-domain partitioned simulation engines

is the 2ð Þ� K þ 1 ð Þ 2K þ 1 composite conversion matrix, known as the Jacobian matrix of the error function F Yð Þ. G and C denote the 2ð Þ� K þ 1 ð Þ 2K þ 1 conversion matrices (Toeplitz) [7]

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 time-domain latency. It must be pointed out that time-domain latency is a phenomenon that is not restricted to heterogeneous circuits. For instance, regimes of operation of pure digital electronic systems typically present a set of variables remaining practically constant within a specific time interval, while others evidence quick variations (fast transitions) in that interval.

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

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–10]. These powerful techniques, denoted as multirate Runge-Kutta (MRK) schemes, split the ODE system of (5) into coupled fast and slow (latent) sub-

, <sup>y</sup>Fð Þ¼ <sup>t</sup><sup>0</sup> <sup>y</sup>F,<sup>0</sup>

, <sup>y</sup>Lð Þ¼ <sup>t</sup><sup>0</sup> <sup>y</sup>L,0,

(19)

or RK schemes) integrate all the DAE, or ODE, unknowns with the same step size.

3.2. Partitioned algorithms for time-step integration

dyFð Þt

dyLð Þt

dt <sup>¼</sup> <sup>f</sup> <sup>F</sup> <sup>t</sup>; <sup>y</sup>F; <sup>y</sup><sup>L</sup>

dt <sup>¼</sup> <sup>f</sup> <sup>L</sup> <sup>t</sup>; <sup>y</sup>F; <sup>y</sup><sup>L</sup>

systems, obtaining

<sup>d</sup><sup>Y</sup> <sup>¼</sup> J Yð Þ¼ G Yð Þþ <sup>j</sup>ΩC Yð Þ (18)

http://dx.doi.org/10.5772/intechopen.75490

29

Numerical Simulation of Electronic Systems Based on Circuit-Block Partitioning Strategies

in which ω<sup>0</sup> ¼ 2π=T is the fundamental frequency. If we substitute (14) into (8), and consider an appropriate harmonic truncation at some order k ¼ K, we will obtain

$$p\left(\sum\_{k=-K}^{+K} Y\_k e^{j k \omega\_0 t}\right) + \frac{d}{dt} \left[q\left(\sum\_{k=-K}^{+K} Y\_k e^{j k \omega\_0 t}\right)\right] = \sum\_{k=-K}^{+K} X\_k e^{j k \omega\_0 t}.\tag{15}$$

The HB method converts the differential problem of (8) into the frequency domain, obtaining the 2ð Þ K þ 1 algebraic equations system

$$\mathbf{F}(\mathbf{Y}) = \mathbf{P}(\mathbf{Y}) + j\Omega \mathbf{Q}(\mathbf{Y}) - \mathbf{X} = \mathbf{0} \tag{16}$$

where Y ¼ ½ � Y�<sup>K</sup>;…; Y0;…; YK <sup>T</sup>, <sup>X</sup> <sup>¼</sup> ½ � <sup>X</sup>�<sup>K</sup>;…; <sup>X</sup>0;…; XK <sup>T</sup>, <sup>j</sup><sup>Ω</sup> <sup>¼</sup> diagð Þ �jKω0;…; <sup>0</sup>;…; jKω<sup>0</sup> .

and P, Q are vectors containing the Fourier coefficients of pyt ð Þ ð Þ and qyt ð Þ ð Þ , respectively. The algebraic equations system of (16) is usually denoted as the harmonic balance system, which can be iteratively solved using the Newton method

$$\mathbf{F}\left(\mathbf{Y}^{[r]}\right) + \frac{d\mathbf{F}(\mathbf{Y})}{d\mathbf{Y}}\Big|\_{\mathbf{Y}=\mathbf{Y}^{[r]}} \left[\mathbf{Y}^{[r+1]} - \mathbf{Y}^{[r]}\right] = \mathbf{0},\tag{17}$$

in which

2.3.2. The harmonic balance method

28 Numerical Simulations in Engineering and Science

will make HB very inefficient.

represented as the Fourier series

x tðÞ¼ <sup>X</sup> þ∞

p X þK

the 2ð Þ K þ 1 algebraic equations system

where Y ¼ ½ � Y�<sup>K</sup>;…; Y0;…; YK

k¼�K

be iteratively solved using the Newton method

F Y½ �<sup>r</sup> � �

þ dF Yð Þ dY

Yke jkω0t !

k¼�∞

an appropriate harmonic truncation at some order k ¼ K, we will obtain

þ d dt <sup>q</sup> <sup>X</sup> þK

Xke jkω0t

The harmonic balance (HB) method [4–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

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 x tð Þ ∈ R satisfying x tðÞ¼ x tð Þ þ T and whose dynamic behavior is described by a unique unknown, y tð Þ∈ R (the generalization to x tð Þ<sup>∈</sup> <sup>R</sup>n, ytð Þ <sup>∈</sup> <sup>R</sup><sup>n</sup> is straightforward). Given that the steady-state regime of the circuit will be periodic with the same period T, both the stimulus and the steady-state solution can be

, ytðÞ¼ <sup>X</sup>

Yke jkω0t

" # !

in which ω<sup>0</sup> ¼ 2π=T is the fundamental frequency. If we substitute (14) into (8), and consider

k¼�K

The HB method converts the differential problem of (8) into the frequency domain, obtaining

and P, Q are vectors containing the Fourier coefficients of pyt ð Þ ð Þ and qyt ð Þ ð Þ , respectively. The algebraic equations system of (16) is usually denoted as the harmonic balance system, which can

<sup>Y</sup>½ � <sup>r</sup>þ<sup>1</sup> � <sup>Y</sup>½ �<sup>r</sup> h i

<sup>T</sup>, <sup>X</sup> <sup>¼</sup> ½ � <sup>X</sup>�<sup>K</sup>;…; <sup>X</sup>0;…; XK

� � � � <sup>Y</sup>¼Y½ �<sup>r</sup> þ∞

Yke jkω0t

<sup>¼</sup> <sup>X</sup> þK

F Yð Þ¼ P Yð Þþ jΩQ Yð Þ� X ¼ 0 (16)

k¼�K

Xke jkω0t

<sup>T</sup>, <sup>j</sup><sup>Ω</sup> <sup>¼</sup> diagð Þ �jKω0;…; <sup>0</sup>;…; jKω<sup>0</sup> .

¼ 0, (17)

, (14)

: (15)

k¼�∞

$$\frac{d\mathbf{F}(\mathbf{Y})}{d\mathbf{Y}} = \mathbf{J}(\mathbf{Y}) = \mathbf{G}(\mathbf{Y}) + j\Omega \mathbf{C}(\mathbf{Y}) \tag{18}$$

is the 2ð Þ� K þ 1 ð Þ 2K þ 1 composite conversion matrix, known as the Jacobian matrix of the error function F Yð Þ. G and C denote the 2ð Þ� K þ 1 ð Þ 2K þ 1 conversion matrices (Toeplitz) [7] corresponding to gyt ð Þ¼ ð Þ dp y t ð Þ ð Þ =dy and cyt ð Þ¼ ð Þ dq y t ð Þ ð Þ =dy.
