IEICE TRANS. FUNDAMENTALS, VOL.E83–A, NO.6 JUNE 2000

1023

PAPER

Special Section of Papers Selected from ITC-CSCC’99

Steady-State Response of Nonlinear Circuits Containing
Parasitic Elements
Takeshi MATSUDA† , Student Member, Yoshifumi NISHIO† , Yoshihiro YAMAGAMI† ,and Akio USHIDA† , Regular Members

SUMMARY We propose here a time-domain shooting algorithm for calculating the steady-state responses of nonlinear RF
circuits containing parasitic elements that is based on both a
modiﬁed Newton and a secant methods. Bipolar transistors and
MOSFETs in ICs have small parasitic capacitors among theirterminals. We can not neglect them because they will gives large
eﬀects to the shooting algorithm at the high frequency. Since our
purpose is to develop a user friendly simulator, we mainly take
into account the relatively large normal capacitors such as coupling and/or by-pass capacitors and so on, because the parasitic
capacitors are usually smaller and contained in the device models. We have developed a very simple simulator only using thefundamental tools of SPICE, which can be applied to relatively
large scale ICs, eﬃciently.
key words: steady-state analysis, RF circuits, time-domain secant method, parasitic capacitors, SPICE

1.

Introduction

It is very important to analyze the steady-state responses for designing communication circuits such as
modulators, mixers, etc. When the attenuation of transient response is suﬃciently large, we can easily calculate the response with the brute-force method. However,
it is sometimes happened that the transient response
continues for a long period due to the small attenuation.
In this case, there are two basic approaches; i.e., thefrequency-domain method [1]–[3] and the time-domain
method [4]–[7]. The former is based on the harmonicbalance method, which is usefully applied to weakly
nonlinear circuits. However, the computational eﬃciency rapidly decreases in the cases when the numberof nonlinear elements increases, and the nonlinearities
become stronger, because the determining equation for
calculating the Fourier coeﬃcients becomes very large
scale. Fortunately, some algorithms for solving large
scale determining equation have been proposed [8], [9]which are eﬃciently applied to the system with thelarge scale sparse Jacobian matrix. Frequency-domain
relaxation method [14] is also eﬃciently applied to theanalysis of relatively large scale weakly nonlinear RF
circuits.
On the other hand, the time-domain method isbased on the transient analysis, where the initial guess
Manuscript received October 4, 1999.
Manuscript revised January 7, 2000.
†
The authors are with the Faculty of Engineering, Tokushima University, Tokushima-shi, 770–8506 Japan.

giving rise to the steady-state response is calculatedby the Newton type shooting method. The ﬁrst one
is based on the Newton Raphson method whose Jacobian matrix is estimated by the analysis of the timevarying sensitivity circuits equal to the number of statevariables [4], [5]. Note that it will take a large computer
memory for a large scale circuit containing many parasitic elements. Fortunately, the transient terms dueto the parasitic elements will be quickly reduced in thesensitivity analysis, so that the Jacobian matrix corresponding to the elements will be approximately replaced by the unit matrix in the Newton iteration [6].
Thus, it is possible to reduce the size of the Jacobianmatrix. Although the extrapolation method is not theoretically guaranteed the convergency, it is a very simple algorithm, and seems to be suitable for the implementation with the SPICE simulator [7].
We propose here an algorithm for calculating thesteady-state responses based on both the modiﬁed Newton and secant methods [10], [11]. Since the secant
method is a type of Newton method whose Jacobianmatrix is successively modiﬁed in the iterations, it canbe applied to calculate both the stable and unstablesteady-state responses if the circuit does not containparasitic elements. Otherwise, our method can be only
applied to caluculate the stable steady-state response.
There have been published many secant methods[19], and we already applied one of them to calculate thesteady-state responses of nonlinear circuits [20] which isbased on a discrete Newton method. The algorithm becomes sometimes unstable near at the solution point. In
this paper, we apply much more stable algorithm using
an orthogonal procedure [10]. Note that although theconvergence ratios of the secant methods are smaller
than the Newton method, they can be usefully applied
to get the solution of a nonlinear system when the Jacobian matrix can not be explicitly obtained.
Thus, the method is suitable for the developmentof a simple SPICE simulator, because we only use thestate variables at every periodic point in the transient
analysis, and need not to use the time-varying sensitivity circuits in the iteration.
At the ﬁrst step in our algorithm, the initial guessis estimated by the modiﬁed Newton method in themeaning that the Jacobian matrix is calculated at thedc operating point, and after then, the matrix is succes-

IEICE TRANS. FUNDAMENTALS, VOL.E83–A, NO.6 JUNE 2000

1024

state response has a period called the total period [3]that is deﬁned by
T =Fig. 1

LSI circuit with the normal capacitors C1 , ..., Cn .

sively modiﬁed by a secant method in the iterations. In
our algorithm, the initial guess giving rise to the steadystate response is calculated for the normal capacitor
voltages after the transient response due to the parasitic capacitors is reduced. The transient period is approximately estimated from the frequency-domain driving point characteristic for the modiﬁed circuit where
the normal capacitors are replaced by the bias voltage
sources. Thus, the characteristic curve can be easily
obtained by the ac-sweep of SPICE, and the poles are
estimated from the approximate rational polynomial
function. We have developed a user friendly simulator using the fundamental tools of SPICE, which canbe applied to the relatively large scale ICs, eﬃciently.
2.

Modiﬁed Newton and Secant Methods

2.1 Eﬀect of Parasitic Elements
To focus on the main idea of our algorithm, consider
a LSI circuit shown in Fig. 1, where C1 , ..., Cn are thenormal capacitors such as coupling and/or by-pass capacitors and so on, which are considered as suﬃciently
large compared to the parasitic capacitors C1p , ..., Cnp ,and they will give large eﬀects to the period of transient response. We call the normal capacitors simply
capacitors, in the following. For example, bipolar transistors have typically 10−12 F order nonlinear capacitors among the emitter-base, the base-collector and theemitter-collector [15], [16]. Depending on circuit conﬁgurations, we may not neglect their eﬀects in the highfrequency such as over 100 MHz. Since MOSFETs also
have small parasitic capacitors such as the gate-drain
CGD , the gate-source CGS , and the bulk-gate, bulkdrain and bulk-source, CBG , CBD , and CBS , respectively [16]. Although they may give an eﬀect to thetransient response in the high frequency region, the effect will be much smaller than that from the normal
capacitors. Furthermore, it is troublesome to consider
them in the circuit simulation, so that we mainly take
into account the normal capacitors in our simulation
algorithm.
Assume that the circuit is driven by multiple inputs e(t) and j(t), whose frequency components are
f1 , ..., fn . If they are the integer relations, the steady-1, ν = 2π/T
G.C.M. {f1 , . . . , fn }(1)

where G.C.M. means the greatest common measure.
Hence, we can estimate the fundamental frequency
component for the multiple frequencies and the period.
Now, let us derive the circuit equation in the formof algebraic-diﬀerential equation. At ﬁrst, let us choose
a normal tree for a given circuit such that it must contain the maximum possible number of the capacitors,
and after then, as many as the parasitic capacitors, andlastly, the resistors. Next, we deﬁne the variables suchas v for the capacitor voltages, v p for the parasitic capacitors and v R for the resistors in the normal tree.
Then, we can describe the circuit equation in thefollowing forms using the fundamental cutset equations
and the some loop equations [17];C(v)v̇ = f 1 (v, v R , v p , νt)(2a)εC p (v p )v̇ p = f 2 (v, v R , v p , νt)(2b)f 3 (v, v R , v p , νt) = 0(2c)where Eqs. (2a) and (2b) are the cutset equations corresponding to the normal tree capacitors and the parasitic capacitors, respectively. Equation (2c) is the loop
equations containing the normal tree resistors and/or
C-E loops [17]. We assume that ε for the parasitic capacitors shows a suﬃciently small constant.
Our time-domain shooting algorithm ﬁnds out thecapacitor voltages v(Tp ) giving rise to the steady-state
response, where Tp is the time such that the eﬀect of theparasitic capacitors becomes negligible in the transient
response. To estimate Tp in a qualitative point of view,
we consider the variational equation at the steady-state
response {v 0 , v R0 , v p0 };v = v 0 + ∆v , v R = v R0 + ∆v R , v p = v p0 + ∆v p (3)Substituting Eq. (3) into Eq. (2), we have


∂f 1
∂f 1
∂f 1


∂vR
∂vp 
C(v 0 )∆v̇ ∂v
f 2 ∂f 2  εC p (v p0 )∆v̇ p  =  ∂∂fv2 ∂∂v

∂vp 
R

0
∂f 3
∂f 3
∂f 3
∂v
∂vR
∂vp


 
∂ C (v )
v̇
∆v
∆v
0
∂v


×  ∆v R  −  ε∂ C p (v p ) v̇ p0 ∆v p 
∂vp
∆v p
0

(4)

where the Jacobian matrix is estimated at the steadystate response. Equation (4) is a time-varying systemwith the period T . From the third row, we have

−1 


∂f 3
∂f 3 ∂f 3∆v∆v R = −
(5)
∆v p
∂v R
∂v ∂v p

MATSUDA et al.: STEADY-STATE RESPONSE OF NONLINEAR CIRCUITS CONTAINING PARASITIC ELEMENTS

1025

Now, assume that the matrices C(v 0 ) and C p (v p0 )are positive and deﬁnite at the steady-state response.
Then, they have the inverse matrices, and substituting
Eq. (5) into Eq. (4), we have the following form;
 



A11 (νt) A12 (νt)
∆v
∆v̇
=
(6)
∆v p
A21 (νt) A22 (νt)ε∆v̇ p
In order to investigate the eﬀect of the parasitic capacitors, we consider the second equation of Eq. (6);ε∆v̇ p = A21 (νt)∆v + A22 (νt)∆v p(7)

First term in the right hand side is a variational value
from the capacitor voltages which behaves like as aforced input in Eq. (7). It is known that the solution ofa linear diﬀerential equation consists of two solutions;
one is the particular solution v pp (t), and another is thegeneral solution v pt (t) satisfying
ε∆v̇ pt = A22 (νt)∆vpt(8)Then, it is known from the Floquet theorem [12] thatthe general solution has a following property;∆v pt (t + T ) = exp(λp T )∆v pt (t),for λp = λ2 /ε (9)where λ2 is a characteristic multiplier of the timevarying diﬀerential Eq. (8) whose real part will be negative. Observe that if |Reλp | is suﬃciently large, thetransient terms due to the parasitic capacitors will bequickly reduced, and the solution of Eq. (7) will quickly
approach to the particular solution v pp (t). It meansthat the response v p (t) only depends on v(t) after thetransient term due to the parasitic capacitors is reduced, and the variable v p (t) no more behaves like as
the state variables. Therefore, after the period Tp , it ispossible to apply our shooting method only to the normal capacitor voltages v(t), and we can ﬁnd the initial
guess giving rise to the steady-state response.
Now, let us estimate the approximate transient period Tp in the case that the cirsuit is at the dc poerating
point. Namely, we calculate the transient response forthe modiﬁed circuit whose the capacitor Cs are replacedby the voltage source v(0) equal to the operating points.
Then, the resultant circuit only contains the parasitic
capacitors, and Tp can be estimate by the transient response. It is also possible to estimate Tp from the frequency response. Namely, let us calculate the driving
point admittance in the frequency domain at any capacitor port† . In practice, the frequency response canbe calculated by the ac-sweep of SPICE at one of thecapacitor ports. The response curve is approximately
described by a rational polynomial as follows;
Y (s) =

a 0 + a 1 s + . . . + a m sm
b0 + b1 s + . . . + bn sn(10)

In order to estimate the approximate transient period
Tp , we need to know the smallest pole in the negativeof the left hand complex plane, so that it is enough toapply the lower order of a rational polynomial function
to the approximation [21].
Remark that if the eﬀect from some of the normal
capacitors to the transient response is negligible, we
can also consider them as parasitic capacitors.
2.2 Our Secant Method
Our purpose is to develop a user friendly simulator using the fundamental tools of SPICE, where we only
consider the capacitor voltages v(t) for determining thesteady-state response. For a large scale circuit, it istroublesome to choose the parasitic capacitor voltages
as the state variables, because the number becomesenormous when the circuit scale increases. Furthermore, they are contained in the device models.
Now, let the period of transient response due to theparasitic elements be Tp . Then, we need to calculate
the initial guesses v(Tp ) giving rise to the steady-state
response, which will satisfy the following determining
equation†† :F (v (Tp )) = v (Tp + T ) − v (Tp ) = 0, for T = 2π/ν(11)
Let us solve the determining Eq. (11) by the iterational
method. It can be eﬃciently calculated by the Newton
method [4] as follows;v

j+1(Tp ) = v (Tp ) −j∂F (v (Tp ))∂v (Tp ) v (Tp )=v j (Tp )× F v j (Tp ) ,j = 0, 1, . . .

−1

(12)

The Jacobian matrix can be calculated by the analysisof time-varying sensitivity circuits for all of the statevariables [4], [5] in the period [Tp , T + Tp ]. Therefore,
the computer eﬃciency will be decreased as the circuit
scale becomes lager and the total period T becomes
longer.
On the other hand, there is a simple algorithm
based on the modiﬁed Newton method, where the Jacobian matrix is estimated at the dc operating point.
Namely, each column of the matrix is calculated by thetransient response to a suﬃciently small impulsive input for corresponding capacitor; ∆v (T ) /∆v (0) for all
the capacitor Cs. If we use the same Jacobian matrix∆v (T )∂F (v(T ))
= I−
(13)
∂v (T )∆v (0)†
It is known that although there are the same numberof diﬀerent driving point admittances as the capacitor Cs,
their poles for all of the admittances are located at the same
points in the complex plane.
††
It is shown in Sect. 2.1 that we can get the exact steadystate response with our secant method after the transientphenomena due to the eﬀect of parasitic capacitors is completely ﬁnished.

IEICE TRANS. FUNDAMENTALS, VOL.E83–A, NO.6 JUNE 20001026

at all of the iterations, the algorithm is called the modiﬁed Newton method. Although the algorithm can beeﬃciently applied to the weakly nonlinear circuits, theconvergence is not guaranteed for the strong nonlinearities.
Hence, we propose here a secant method [10] suchthat the Jacobian matrix [11] is successively improvedin the latter iterations. Let J 0 be the initial Jacobianmatrix estimated by the above method. Then, we haveat the ﬁrst iteration
J 0 δv 0 (Tp ) = −F v 0 (Tp )
1

0

(14)
0

for v (Tp ) = v (Tp ) + δv (Tp )Thus, we need to estimate the modiﬁed Jacobianmatrix J 1 . Suppose it satisﬁes the following relation;
F v 1 (Tp ) = F v 0 (Tp ) + J 1 δv 0 (Tp )(15)where J 1 = J 0 + D 0
D 0 corresponds to the variational Jacobian matrix to
be determined. From Eqs. (14) and (15), we have010

F v (Tp ) = D δv (Tp )(16)

where D 0 can be solved as follows;
D0 =

x

(17)

T
(x0 ) δv 0 (Tp )T

where x0 is an arbitrary vector satisfying x0 δv 0 (Tp )= 0. In the same manner, we have the following iteration;
δv j (Tp ) = − J j−1

F (v j (Tp ))v j+1 (Tp ) = v j (Tp ) + δv j (Tp )j

D =F v j+1 (Tp )xj

J j+1 δv j (Tp ) = F (v j+1 (Tp )) − F (v j (Tp ))
for J

(19)

(20)

T

j

(18)

T

(xj ) δv j (Tp )

j+1

(21)

j

= J + D , j = 1, 2, . . .

where if j > n, xj is chosen orthogonal to the previous
n−1 steps [18], δv j−n+1 (Tp ), . . . , δv j−1 (Tp ), and if j <n, we can only demand that xj is orthogonal to theavailable j − 1 steps, δv 1 (Tp ), . . . , δv j−1 (Tp ). Observethat our secant method only uses the data v(t) at t =Tp and t = T + Tp in the transient response, and it isobtained with the SPICE. Thus, we can easily develop auser friendly simulator only using fundamental tools of
SPICE, where we need not to derive any troublesome
circuit equations. It makes our scant method much
more powerful.
Our secant algorithm

S.2 Calculate the solution at j + 1st iteration
v j+1 (Tp ) = v j (Tp ) + δv j (Tp )
S.3 If F v j+1 (Tp )  < ε for a given small ε, thenstop. Otherwise, calculate the variational value of
the Jacobian matrix
Dj =F v j+1 (Tp )

xj

T

T

(xj ) δv j (Tp )S.4 Set J j+1 = J j + D j .
S.5 Set j = j + 1 and go to Step 1.

0 T

F v 1 (Tp )S.0 Set j = 0, and set the initial guess v 0 (Tp ) equal tothe bias voltages. The initial Jacobian matrix J 0
is estimated by Eq. (13) at the dc bias points.
Applying the rational approximation to the driving point admittance obtained by the ac-sweep of
SPICE, we estimate the approximate transient period Tp due to the parasitic elements. Set x0 =[1, 0, . . . , 0], and stopping condition ε.
S.1 Calculate the transient response in the period of[0, T + Tp ] from the initial guess v(0) = v j (Tp ).
Let us estimate the variation by −1
δv j (Tp ) = − J j
F v j (Tp )In this algorithm, the vector xj must be chosen orthogonal to the previous vectors δv i , i = j − n + 1, . . . , j − 1,which is eﬃciently executed with the Schmidt orthonormalization procedure [18]. Observe that, to implement
our algorithm, we only need to execute the transient
analysis in the period [0, Tp + T ] and get v(Tp ) andv(T + Tp ).
The ﬂow chart of our secant method is shown in
Fig. 2. In our algorithm, the data v j (Tp ) and v j (Tp +T ) are obtained from the transient analysis of SPICE
and they are transferred to C-language program, where
the initial guess v j+1 (Tp ) is calculated with the secant
method. After then, the data from C-language programis again transferred to SPICE, and so on.
Remark 1: It is known that if one of the variables inδv j (Tp ) = v j+1 (Tp )−v j (Tp ) approaches to zero, the xj
obtained by the orthonormalization procedure [18] may
have serious error. In this case, the iteration sometimes
happens to become unstable. Thus, we will recommend
to remove the variable vi from v if|δvij | < δ, for some ifor a suﬃciently small δ, where δ depends on both thedigit of computer and the truncation error in the numerical integration method [13]. Thus, the dimensionof J j and D j should be reduced by one, and the orthogonal vector xj should be estimated for the remaining
variables.
Remark 2: If we deﬁne the convergence ratio by p in
the following relation

MATSUDA et al.: STEADY-STATE RESPONSE OF NONLINEAR CIRCUITS CONTAINING PARASITIC ELEMENTS

1027

(a)

Fig. 2||v j+1 (Tp ) − v̂(Tp )|| = k||v j (Tp ) − v̂(Tp )||p
where v̂ means the exact solution, then, the convergence rate of the secant method is p = 1.62 [11]. Although the ratio is smaller than the Newton method
p = 2, our secant algorithm is a very simple because
we need not solve the sensitivity circuit. Therefore,
our method seems to be suitable for developing a user
friendly simulator using SPICE.
3.

(b)

Flow chart of the secant method.

Illstrative Examples

3.1 RLC Circuit with Nonlinear ResistanceIn order to investigate the eﬀect of parasitic capacitors
in the steady-state solution, we consider a simple RLC
circuit as shown in Fig. 3(a), where the parasitic capacitor Cp is chosen 10% of the normal capacitor C. Since
λp = −20 due to the parasitic capacitor is very large
in negative, the eﬀect in the transient response will be
negligible after 1 [sec].
Although the number of state variables is 3 in thisexample, the transient phenomena will behave like as
2 state-variables circuit after the transient period due
to Cp is reduced. Thus, we can ﬁnd out the initial
guess v(Tp ) giving rise to the steady-state solution after
Tp = 1 [sec].
In Table 1 where the exact solution is obtainedby the Newton method [4]. The second result is obtained by solving the determining equation F (v(0)) =Fig. 3 (a) RLC circuit with nonlinear resistance. R = 1, L =3 e(t) = sin 2πt. (b)1, C = 1, Cp = 0.1, iR = vR + vR
Convergence rate. Error = F (v j ) .

Table 1
Exact solution
F(v(0)) = 0
F(v(Tp )) = 0

Comparisons of the initial guesses.
vC (Tp )−7.7357 × 10−3
−5.8436 × 10−3
−7.7355 × 10−3

iL (Tp )−0.16130−0.16196−0.16130vCp (Tp )1.4537-v(0) − v(T ). In this case, even if the convergence algorithm has a solution, it will be the fault solution inthe (vC (0), iL (0), vCp (0))-plane† . The third result is obtained with our secant method, where we have chosen
Tp = T for simplicity. After the transient due to Cp is
reduced, the circuit will behave like as in the 2 dimensional plane. Therefore, the result is almost equal tothe exact solution. The convergence ratio is suﬃcientlarge as shown in Fig. 3(b).
3.2 RC Ampliﬁer
Now, consider a simple RC ampliﬁer circuit shown in
Fig. 4(a). The transistor has the small parasitic capacitors among the emitter, base and collector [15], which†
Note that there will exist an inﬁnite number of solutions satisfying F (v (0)) = v (0) − v (T ) in 3 dimensional (vC (0), iL (0), vCp (0))-plane because the algorithm using F (v (0)) is in 2-dimensional plane.

IEICE TRANS. FUNDAMENTALS, VOL.E83–A, NO.6 JUNE 2000

1028

F (v(0)) neglecting the eﬀect of the parasitic capacitors never converges to any solution. Note that ouralgorithm using F (v(Tp )) can ﬁnd out the steady-statesolution within 28 [sec], and the transient analysis with
the SPICE gets the solution within 241 [sec].
3.3 Four Phase Mixer(a)This is an example of a relatively large scale four phasemixer circuit that consists of 122 bipolar transistors,
and some capacitors, as shown in Fig. 5(a). In this case,the time-domain shooting method based on the Newton
method [4] will be very time-consuming, because thenumber of the parasitic is more than 250 in this circuitand the total period for two inputs is large.
In designing of the mixer circuit, it is very important to know the intermodulated frequency components
in the output waveform. Assume that the circuit isdriven by two signals as follows;
v1 (t) = 282 cos(2π × 13 × 106 t) [mV] : Local oscillator(b)Fig. 4 (a) RC ampliﬁer. R1 = 490 kΩ, R2 = 2 kΩ, R3 =100 kΩ, R4 = 200 Ω, R5 = 600 Ω, C1 = 10 nF, C2 =10 nF, C3 = 10 nF, e(t) = 0.1 sin 2π108 t [V]. (b) The driving
point characteristic curve at port C1 .

cannot be neglected at higher than 100 [MHz].
At ﬁrst, we choose the coupling capacitors C1 , C2
and the by-pass capacitor C3 as the normal capacitors.
To investigate the eﬀect of the parasitic capacitors, thenormal capacitors are replaced with the operating voltage sources v1 (0) = 1.007 [V], v2 (0) = 6.641 [V] andv3 (0) = 0.338 [V], which are obtained by the dc analysis of SPICE. Now, we apply the ac-sweep of SPICE,
and get the driving point characteristic curve as shownin Fig. 4(b). We approximate it with the second orderrational function, because we enough to know only thesmallest pole in negative [21]; i.e.,
Y (s) =

3.27×10−4 + 2.81×10−11 s + 9.48×10−20 s2
1. + 1.415 × 10−9 s

It is well-known from the circuit theory that althoughfrequency responses from the other ports have the different characteristics, their poles are located in almost
same point;
λp = −7.1 × 108
Therefore, we can hope to get the exact solution with
Tp = 1./7.1 × 10−8  10−8 [sec], andF (v(Tp )) = v(Tp ) − v(2Tp ), for Tp = T
Remark that we can get the solution in 7 iterations.
On the other hand, the time-domain secant method

v2 (t) = 9.12 sin(2π × 14 × 106 t) [mV] : Input signal
The fundamental frequency of the mixer outputs
are 1 [MHz] which are obtained at the v1out , ..., v4out
terminals in Fig. 5(a). In this case, C1 , C2 and C3 arecoupling capacitors between the sub-circuits which will
give large eﬀect on the transient phenomena. On theother hand, the capacitor Cs used as ﬁlter circuits inthe output stage do not give the eﬀect to the righthand sub-circuits, because they are separated by the
buﬀer ampliﬁers. Therefore, we have chosen the voltages of C1 , C2 and C3 as the state variables in our
secant method.
Furthermore, since the lower side sub-circuits inthis circuit are also separated by a buﬀer ampliﬁer fromthe upper one, we can assume that the voltages at C1 ,C2 and C3 will only contain the fundamental frequencyof ω1 = 2π × 13 × 106 and its higher harmonics. This
makes the analysis much easier, because we can deﬁne
the total period T = 1/(13 × 106 ) [sec].
Now, let us apply our algorithm to the circuit.
Since the transient period due to the parasitic capacitors is very short compared to the input frequency 13 [MHz], we can get the same result for both
F (v(0)) = 0 and F (v(T )) = 0 in our algorithm. The
convergence rate for F (v(T )) = 0 is shown in Fig. 5(b).
We found from the result that the convergence ratios inthe ﬁrst 3 iterations are smaller compared to the following iterations. The result can be explained as follows
that the Jacobian matrix has 3 × 3 elements and theinitial matrix J 0 estimated at the dc operating point
does not seem to be a good approximation. However,the Jacobian is successively improved in the iterations.
Thus, the convergence ratio after the ﬁrst 3 iteration
becomes larger.

MATSUDA et al.: STEADY-STATE RESPONSE OF NONLINEAR CIRCUITS CONTAINING PARASITIC ELEMENTS

1029

(a)

(c)

(b)

(d)
Fig. 5 (a) Four phase mixer circuit. C1 = C2 = 20 pF, C3 = 2.5 nF, C = 1 µF. (b)Convergence rate. Error = F (v j ) . (c) Steady-state waveform with our method. (d)Transient waveform obtained from SPICE.

We compared the computational times for otherdiﬀerent methods.
In Table 2 where “PSS” means the periodic steadystate analysis based on the time-domain shooting
method, whose SpectreRF is widely used in the RF circuit simulation. We can conclude that the SpectreRF(PSS) in this example is ineﬃcient compared with our
method. Thus, we found that our secant method can beapplied to the steady-state analysis of relatively large
scale ICs, eﬃciently. The output waveforms obtained
from our method and transient analysis of SPICE are
shown in Figs. 5(c), (d), respectively.

IEICE TRANS. FUNDAMENTALS, VOL.E83–A, NO.6 JUNE 2000

1030
Table 2Comparison of the computational times.

PSPICE (Transient)
SpectreRF (PSS)
Relaxation method [14]
Our secant method

Computation time704.21 [sec]140.54 [sec]54.99 [sec]24.48 [sec]Remark that the total period of the circuit is deﬁned by the diﬀerence between the local oscillator andthe input signal frequencies. Thus, when the diﬀerence is very small, the period becomes very longer,
and the time-domain shooting method using SpectreRF becomes time-consuming. Note that we can not
get the steady-state response for f1 = 13 [MHz] andf2 = 13.1 [MHz] with our 160 Mbytes computer becauseof the memory over [14]. The computational eﬃciencyof our method for this example does not be changed
even if the total period becomes longer.
4.

Conclusions and Remarks

In this paper, we have shown an eﬃcient time-domainsecant method for calculating the steady-state response. Although the convergence rate is smaller thanthe Newton method, the algorithm is very simple andsuitable for the development of the user friendly simulator. In our simulator, the initial guess from the
SPICE is improved by the secant method written by
the C-language program, and the initial guess is again
returned to the SPICE, and so on. We continue thesame iteration until the exact solution can be obtained.
Thus, our algorithm is very simple, and need not derive
any troublesome circuit equations and the sensitivity
circuit analysis. Our simulator will be eﬃciently applied to relatively large scale RF circuits such as modulators and mixers. In this paper, we assumed that theeﬀect of parasitic capacitors is smaller than the normal capacitors. However, the eﬃciency of our method
may be decrease if the transient response Tp due to theparasitic elements becomes longer.
Acknowledgement
The authors would like to express their gratitude to
Mr. Kimihiro Ogawa and Mr. Masayuki Takahasi atSony Co. for useful discussions and suggestions concerning the algorithm and results.
References[1] K. Kundert, J. White, and A. Sangiovanni-Vincentelli,
Steady-State Methods for Simulating Analog and Microwave Circuits, Kluwer Academic Publishers, Boston,
1990.
[2] H.G. Brachtendorf, G. Welsch, and R. Laur, “Fast simulation of the steady-state of circuits by the harmonic balancetechnique,” Proc. IEEE Int. Symp. on Circuits & Systems,
pp.1388–1391, 1995.

[3] A. Ushida, T. Adachi, and L.O. Chua, “Steady-state analysis nonlinear circuits based on hybrid method,” IEEE
Trans. Circuits & Syst.-I: Fundamental Theory and Applications, vol.39, no.8, pp.649–661, 1992.
[4] T.J. Aprille, Jr. and T.N. Trick, “Steady-state analysis ofnonlinear cirsuits with periodic input,” Proc. IEEE, vol.60,
pp.108–114, 1972.
[5] L.O. Chua and A. Ushida, “Algorithm for computing almost periodic steady-state response of nonlinear systems tomultiple input frequencies,” IEEE Trans. Circuits & Syst.,
vol.CAS-28, pp.953–971, 1981.
[6] K. Kakizaki and T. Sugawara, “A modiﬁed Newton method
for the steady-state analysis,” IEEE Trans. Comput.-Aided
Des. Intergrated Circuits & Syst., vol.CAD-4, pp.354–360,
1985.
[7] S. Skelboe, “Computation of the periodic steady-state response of nonlinear networks by extrapolation methods,”
IEEE Trans. Circuits & Syst., vol.CAS-27, pp.161–175,
1980.
[8] R. Telichevesky, K. Kundert, and J. White, “Eﬃcient
steady-state analysis based on matrix-tree Krylovsubspacemethod,” Proc. Int. Conf. on Computer-Aided Design, San
Francisco, DAC’95, pp.480–484, 1995.
[9] R. Telichevesky, K. Kundert, and J. White, “Eﬃcient AC
and noise analysis of two-tone RF circuits,” Proc. Int. Conf.
on Computer-Aided-Design, Las Vegas, pp.292–297, 1996.
[10] J.G.P. Barnes, “An algorithm for solving nonlinear equation
based on the secant method,” Computer Journal, pp.66–72,
1965.
[11] J.F. Traub, Iterative methods for the solution of equations,
Prentice-Hall, Englewood Cliﬀs, N.J., 1964.
[12] J.K. Hale, Oscillations in Nonlinear Systems, McGraw-Hill,
1963.
[13] L.O. Chua and P.-M. Lin, Computer Aided Analysis ofElectronic Circuits: Algorithms and Computational Techniques, Prentice-Hall, Englewood Cliﬀs, N.J., 1975.
[14] Y. Yamagami, Y. Nishio, A. Ushida, M. Takahashi, andK. Ogawa, “Analysis of communication circuits based on
multidimensional Fourier transformation,” IEEE Trans.
Comput.-Aided Des. Intergrated Circuits & Syst., vol.18,
pp.1165–1177, 1999.
[15] H.K. Gummel and H.C. Poon, “An integral charge control
model of bipolar transistors,” Bell Syst. Tech. J., vol.94,
pp.827–852, 1970.
[16] R.L. Geiger, P.E. Allen, and N.R. Strader, VLSI Design
Techniques for Analog and Digital Circuits, McGraw-Hill,
New York, 1990.
[17] R.A. Rohhrer, Circuit Theory: An Introduction to theState Variable Approach, McGraw-Hill, 1970.
[18] L.A. Zadeh and C.A. Desoer, Linear System Theory: The
State Space Approach, McGraw-Hill, 1963.
[19] J.M. Ortega and W.C. Rheinboldt, Iterative Solution of
Nonlinear Equations in Several Variables, Academic Press,1970.
[20] A. Ushida, “A steady-state analysis of nonlinear dynamic
systems,” IECE Trans., vol.J61-A, no.3, pp.231–238, 1978.
[21] E. Chiprout and M.S. Nakhla, “Analysis of interconnect
networks using complex frequency hopping (CFH),” IEEE
Trans. Comput.-Aided Des. Integrated Circuits & Syst.,
vol.14, no.2, pp.186–200, 1995.

MATSUDA et al.: STEADY-STATE RESPONSE OF NONLINEAR CIRCUITS CONTAINING PARASITIC ELEMENTS

1031
Takeshi Matsuda
was born in Wakayama, Japan, on 1976. He recieved the
B.E. degrees from Tokushima University,
Tokushima, Japan, in 1999. He is currently working towards M.E. degree at thesame university. His research interest is incomputer-aided analysis of nonlinear systems.

Yoshifumi Nishio
received the B.E.,
M.E., and Ph.D. degrees in ElectricalEngineering from Keio University, Yokohama, Japan, in 1988, 1990 and 1993, respectively. In 1993, he joined the Department of Electrical and Electronic Engineering at Tokushima University, Tokushima, Japan, where he is currently an
Associate Professor. His research interests are in chaos and synchronization phenomena in nonlinear circuits. Dr. Nishio
is a member of the IEEE.

Yoshihiro Yamagami
graduatedfrom the Department of Electrical Engineering, Anan Technical College, Tokushima, Japan, in 1982. He joined the Department of Electrical and Electronic Engineering of Tokushima University, Tokushima, Japan, in 1984. His research interests are computer-aided simulations of
nonlinear circuits and high-speed interconnect cirsuits.

Akio Ushidareceived the B.E.
and M.E. degrees in electrical engineering from Tokushima University in 1961
and 1966, respectively, and the Ph.D. degree in electrical engineering from University of Osaka Prefecture in 1974. He wasan associate professor from 1973 to 1980
at Tokushima University. Since 1980 hehas been a Professor in the Department
of Electrical Engineering at the university.
From 1974 to 1975 he spent one year as
a visiting scholar at the Department of Electrical Engineeringand Computer Sciences at the University of California, Berkeley. His current research interests include numerical methods
and computer-aided analysis of nonlinear systems. Dr. Ushida isa member of IEEE.