Extended Summary

本文は pp.543-550

Spice-Oriented Algorithm for Peak Search and Stability Assessment
for Frequency Response
Hiroshige Kataoka Student Member
Yoshihiro Yamagami Non-member
Yoshifumi Nishio Non-member
Akio Ushida Non-member(Dept. Electrical and Electronic Eng.,
Tokushima University,
hiroshige@ee.tokushima-u.ac.jp)(Dept. Electrical and Electronic Eng.,
Tokushima University)(Dept. Electrical and Electronic Eng.,
Tokushima University,
nishio@ee.tokushima-u.ac.jp)(Dept. Electrical and Electronic Eng.,
Tokushima University)

Keywords: Spice, nonlinear circuit, harmonic balance method, peak search, stability assessment
In this study, we propose a Spice-oriented algorithm forpeak search and stability assessment for the frequency response. We obtain the frequency response by using the sinecosine circuit, Fourier transformation circuit and solution
tracing-curve circuit, which have been developed in our previous studies. We combine the idea using a diﬀerentiator anda nonlinear limiter with them in order to search the exact
peaks and propose a stability assessment algorithm based on
the Floquet theory. In our algorithm, we can simulate almost
all tasks with Spice.
In this article, we introduce our algorithm to the circuitthat is shown in Fig. 1.
Figure 2 shows the simulation result of the frequency response obtained by the proposed algorithm. our algorithm
is very eﬀective to search exact peak values of the frequency
response curve. The obtained peaks are summarized in Table 1.

We assess the stability of the solutions at the points thatare marked in Fig. 2. Table 2 shows the calculated eigenvalues of Φ(T ) for the three cases which obtained by using ouralgorithm.
We can see that the point (ω, i1 ) = (4.570, 1.902) is stable,
because all of eigenvalues satisfy |λ| < 1. However, for theother two points, (ω, i1 ) = (4.927, 0.628) and (5.560, 2.000),
the solutions are unstable, because some eigenvalues do notsatisfy |λ| < 1.
The simulation results for example circuits showed thatthe proposed algorithm worked well.

Fig. 2. Frequency response obtained by the proposed algorithm.
Fig. 1.

Table 2.

Example circuit with strong nonlinearity.

Table 1.

λ1
λ2
λ3
λ4
λ5

Frequencies and voltages at peaks.

ω[rad/s]
i1 [A]

peak1
5.067
2.4545peak2
6.151
2.5000peak3
6.420
3.4450–8–

ω = 4.570−0.622 + j0.554−0.622 − j0.554
0.495 + j0.625
0.495 − j0.625
0.8394

Eigenvalues of Φ(T ).
ω = 4.927
0.726 + j0.825
0.726 − j0.825
0.373 + j0.499
0.373 − j0.499
0.8501

ω = 5.560
2.530
0.545 + j0.638
0.545 − j0.638
0.86620.2862

PaperSpice-Oriented Algorithm for Peak Search and Stability Assessment
for Frequency Response
Hiroshige Kataoka∗
Yoshihiro Yamagami∗
Yoshifumi Nishio∗
Akio Ushida∗

Student Member
Non-member
Non-member
Non-member

In this study, we propose a Spice-oriented algorithm for peak search and stability assessment for the frequency response. We obtain the frequency response by using the sine-cosine circuit, Fourier transformation
circuit and solution tracing-curve circuit, which have been developed in our previous studies. We combine theidea using a diﬀerentiator and a nonlinear limiter with them in order to search the exact peaks and proposea stability assessment algorithm based on the Floquet theory. In our algorithm, we can simulate almost all
tasks with Spice.
Keywords: Spice, nonlinear circuit, harmonic balance method, peak search, stability assessment

ing equations of the HB method by using Spice. In thissubsection, we explain how we can derive the sine-cosine
circuit for simple passive element cases.
First, we set a voltage and a current with Fourier series;⎧n⎪⎪ v = V0 +(Vsk sin kωt + Vck cos kωt)
⎪
⎪
⎨
k=1
· · · (1)n⎪⎪⎪⎪ i = I0 +(Isk sin kωt + Ick cos kωt)⎩

1. Introduction
Spice is a very convenient tool for circuit simulationand is used by many researches. However, when weanalyze nonlinear circuits, we have to make consideration for harmonic components. Because the frequency
response is distorted by the inﬂuence of the nonlinear
elements. In that case, the Harmonic Balance (HB)
method (1)–(3) is useful and we have developed Spiceoriented HB method in our past studies (4)–(6) .
When we simulate nonlinear circuits, it is very important to ﬁnd out the locations and the frequencies
giving large peak values of the voltages. If the circuit
includes some nonlinear elements, resonance curve becomes complicated. In that case, the stability assessment is also very important. In this article, we proposea Spice-oriented algorithm for peak search and stability
assessment for frequency response of nonlinear circuits.
Section 2 reviews the Spice-oriented HB method algorithm for nonlinear circuits. Section 3 shows the peaksearch algorithm for frequency response. Section 4 shows
the algorithm for stability assessment based on the Floquet theory (3) (7) . Illustrated examples of the proposed
algorithm for the peak search and the stability assessment are shown in Sec. 5. Section 6 gives discussions
and Sec. 7 concludes the article.

k=1

A current through a capacitor is given by
i=C

From Eqs. (1) and (2), we can express the current asi=n(−kωCVck sin kωt + kωCVsk cos kωt).

k=1· · · · · · · · · · · · · · · · · · · · · (3)From Eq. (3), we can express the relation between thecoeﬃcients of sine and cosine components as follows;
Isk = −kωCVck· · · · · · · · · · · · · · · · · · · · · · · · · (4)Ick = kωCVsk

2. Spice-Oriented HB Method

In the case of an inductor, we can express the voltage
across an inductor as

In this section, we review the algorithm of Spiceoriented HB method, which has been developed by theauthors (4)–(6) .
2.1 Sine-Cosine Circuit
The sine-cosine circuit has been introduced in order to solve the determin-

v=

n


(−kωLIck sin kωt+kωLIsk cos kωt) · · · (5)k=1

where the current through an inductor is given by∗ Dept. Electrical and Electronic Eng.,

Tokushima University,
2-1, Minami-Josanjima, Tokushima 770-8506
c 2011 The Institute of Electrical Engineers of Japan.
dv. · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · (2)
dt

v=L
543

di
. · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · (6)dtFig. 1.

Example RLC circuit.

Fig. 3.

Fig. 2.

Fourier transformation circuit.

⎧
 2π
1
⎪
⎪
V
=
f (q)d(ωt)⎪
0⎪⎪
2π  0
⎪
⎪
2π
⎪
1
⎨
Vck =
(f (q) cos kωt)d(ωt)· · · · · · · · (10)
π
0

⎪
⎪
⎪
1 2π
⎪
⎪
⎪
⎪ Vsk = π 0 (f (q) sin kωt)d(ωt)⎪⎩(k = 1, 2, . . . , n)Sine-cosine circuit for Fig. 1

Equations for coeﬃcient of sin kωt and cos kωt are givenby
Vsk = −kωLIck· · · · · · · · · · · · · · · · · · · · · · · · · (7)Vck = kωLIsk

In Spice, the coeﬃcients cannot be given in the analytical forms. Therefore, we apply the trapezoidal formula.
 bhg(ωt)d(ωt) = (g(ωt0 ) + g(ωt2M ))2a+ h(g(ωt1 ) + g(ωt2 ) + · · · + g(ωt2M−1 )),ωtm = 0, π/M, ..., (2M − 1)π/M, 2π, · · · · · (11)As an example, we consider a simple linear RLC circuit as Fig. 1. In this circuit, coeﬃcients of sin kωt andcos kωt for the capacitor and the inductor are given by
Eqs. (4) and (7), respectively. We make the sine-cosine
circuits satisfying Eqs. (4) and (7) to obtain the valuesof the coeﬃcients as shown in Fig. 2.
In Fig. 2, the capacitor is replaced by coupled voltagecontrolled current sources and the inductor is replacedby coupled current-controlled voltage sources in the sinecosine circuit.
2.2 Fourier Transformation Circuit
When
we analyze nonlinear circuits with Spice-oriented HB
method, we introduce the Fourier transformation circuit in order to obtain the nonlinear characteristics ofthe elements including many harmonic components.
Let us consider the nonlinear capacitor whose characteristics are given by

where g(ωt) is the function of (ωt), h = (b − a)/2M , andM is the number of the division for an interval (0, π).
From Eqs. (8), (9), (10) and (11), we can obtain theequations for the coeﬃcient of sin kωt and cos kωt as

1 2π
Vck =
(f (q) cos kωt)d(ωt)
π 0
11(f (q0 )+f (q2M ))+ (f (q1 ) cos kωt1=
2MM
+f (q2 ) cos kωt2 +. . .+f (q2M−1 ) cos kωt2M−1 ), 2π· · · · · · · · · · · · · · · · · · · · · (12)1(f (q) sin kωt)d(ωt)π 0
1(f (q1 ) sin kωt1 + f (q2 ) sin kωt2=
M
+ · · · + f (q2M−1 ) sin kωt2M−1 ) · · · · · · · · (13)Vsk =v = f (q). · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · (8)Suppose the input and the output waveforms are written
as follows;
⎧
n

⎪
⎪
⎪ q = Q0 +(Qsk sin kωt + Qck cos kωt)⎪⎨
k=1
· · · (9)n⎪⎪⎪(Vsk sin kωt + Vck cos kωt)⎪
⎩ v = V0 +where f (qm ) = f (qωtm ).
Fig. 3 shows the Fourier transformation circuit. The
circuit is realized to satisfy Eqs. (12) and (13).
In the Fourier transformation circuit, we obtain the
Fourier expansion as voltage values by describing thenonlinear characteristics with ABM (Analog Behavior
Model). Namely, if we can describe the characteristics with ABM, the nonlinear model can be handled in

k=1

The coeﬃcients can be calculated by the Fourier series
expansion as follows;
544

IEEJ Trans. EIS, Vol.131, No.3, 2011

Spice-Oriented Algorithm for Peak Search and Stability Assessment

Fig. 4.

Solution tracing circuit.
Fig. 5.

the proposed method even for complicated characteristics. However, it may be better to develop a software toderive the Fourier transformation circuit automatically,
if we often handle such complicated models. Further,
if we perform a multi-tone analysis, we have to introduce two dimensional Fourier transformation circuit. It
is theoretically possible with a similar concept to the
one-dimensional case, although the circuit becomes verylarge in order to treat a huge number of frequency components.
2.3 Solution Tracing Circuit
Even we combine the sine-cosine circuit and the Fourier transformation circuit, we cannot obtain unstable region of the frequency curve, because we set the frequency as time in
Spice to solve the sine-cosine circuits. In this section,
we explain the Solution Tracing Circuit (STC) (8) (9) realizing the arc-length method (10) (11) .
First, we express the ark length in (n + 1) dimensional
euclidean space as Eq. (14)ds =Fig. 6.

on the frequency characteristics curve. Hence, v(ω)
needs to be ﬁrstly diﬀerentiated. In order to ﬁnd theexact peak points, the output of the diﬀerentiator islimited and expanded with a nonlinear limiter, which
consists of a limiter and diodes as shown in Fig. 5.
We suppose the characteristics of the limiter as follows;⎧⎨ −Vmax : for vin < −VL
va = K vin : for − VL ≤ vin ≤ VL .· · · (17)⎩Vmax : for VL < vin(dx1 )2 + (dx2 )2 + · · · + (dxn+1 )2 · · · · (14)

i=1

dvi
dt

2

+

dvω
dt

2

= 1 · · · · · · · · · · · · · · · · · (15)In our example, we set VL = 0.1 [mV ], Vmax = 10 [V ],and K = 1.0 × 105 . Also, the current through the nonlinear resistor of the diodes is assumed to be given byIs exp(λva ): for va > 0
iout =· · · · · (18)−Is exp(−λva ) : for va < 0

where vi (i = 1, 2, . . . , n) are coeﬃcients for the valuesof the voltages and the currents in Eqs. (4) and (7) andvω corresponds to ω. They are realized by using diﬀerentiators (simply realized by capacitors with 1 (F ) in
Spice). The circuit in Fig. 4 realizes to satisfy the arclength method Eq. (15). In this circuit, the voltages corresponding to the coeﬃcients are inputted after diﬀerentiated and squared via the voltage-controlled current
source (VCCS). If we set the voltage of the node a asv̇ω , Iω = v̇ω2 can be obtained by multiplier and voltagecontrolled current source VCCS (MVCCS). Then, thenode voltage v̇ω is integrated to obtain vω . Note that R
in Fig. 4 is a very large resistor used only to avoid the
L − J cut-set.

where Is = 10−12 and λ = 40. Then, the output characteristics of the nonlinear limiter is obtained as shownin Fig. 6.
The expansion factor of K is set to be large enough.
By virtue of the eﬀect of the nonlinear limiter, the transient analysis around the zero points of the input is executed with a very small step size in Spice. As a result,
we can obtain the exact peak values of the curve automatically.
4. Assessment of Stability

3. Spice-Oriented Peak Search Algorithm

In this section, we explain the proposed peak stability
assessment algorithm for frequency response.
We suppose that the circuit equation is given as

In this section, we explain the proposed peak searchalgorithm for frequency response.
In order to search the exact peaks of the frequency
characteristics curve obtained by using the abovementioned algorithm, we use a diﬀerentiator and a nonlinear limiter. The frequencies ω at the peak voltages
satisfy
電学論 C，131 巻 3 号，2011 年

Characteristics of nonlinear limiter.

dv(ω)= 0, · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · (16)dωIn order to trace the solutions curve by Spice, we replace the diﬀerentiation with respect to the arc-length sby the time t. Also we replace xi by voltages vi in Spice.
From this, we can obtain Eq. (15).
nNonlinear limiter.

f (ẋ, x, y, ωt) = 0, · · · · · · · · · · · · · · · · · · · · · · · · (19)and make the variational equation for the regular periodic solution of (x̂, ŷ). First, we assume the small change
545

quantity as (Δx, Δy) as
x = x̂ + Δxy = ŷ + Δy

· · · · · · · · · · · · · · · · · · · · · · · · · · · · (20)and substitute Eq. (20) to Eq. (19). We obtain the equation as

⎤
⎡
˙

Δx

∂f
∂f
∂f
˙ x̂, ŷ, ωt) +⎦⎣f (x̂,= x̂ Δx = 0.
∂ ẋ ∂x ∂y  x
y = ŷ
Δy· · · · · · · · · · · · · · · · · · · · · (21)Fig. 7. LRC ladder circuit including nonlinear
capacitors.

In Eq. (21), the ﬁrst term is regular periodic solution
and the second term is variational equation. We change
the second term as
˙ = A(t)Δx. · · · · · · · · · · · · · · · · · · · · · · · · · · · · · (22)Δx
In Eq. (22), A(t) is the periodic function. We apply
the Floquet theory for this periodic function. We write
the elementary matrix of the periodic solution as Φ(t).
From this reason, the solution after one period from initial value of Δx(0) is given as follows;

Fig. 8.

Sine-cosine circuit for Fig. 7.

Δx(T ) = Φ(T )Δx(0). · · · · · · · · · · · · · · · · · · · · · · (23)Hence, when the eigenvalues λk (k = 1, 2, . . . , n) of Φ(T )satisfy |λk | < 1 for all k, the regular periodic solution x̂
can be assessed as stable.
In this study, we derive the variational circuit whichcorresponds to the variational equation and perform thetransient analysis of Spice just for one period in orderto obtain the components of the elementary matrix. We
should repeat the transient analysis by giving diﬀerent
initial conditions to obtain all the components numerically. However, the number of the repeat is at most thesame as the number of the state variables of the circuit.
Further, it should be mentioned that we do not have tochange the structure of the variational circuit even whenthe voltages of the regular periodic solution are changed.

Fig. 9. Frequency response obtained by the conventional method in Sec. 2 (circuit in Fig. 7).

5. Illustrative Examples
5.1 Example 1
As the ﬁrst example, we consider
an LRC ladder circuit including nonlinear capacitors as
Fig. 7. We consider the case that all inductors are linear
and have the same inductance as L1 = L2 = L3 = L4 =
0.1 [H], and that all nonlinear capacitors have the same
nonlinear characteristics as G(q) = q + 0.8q 3 . The otherparameters are set as R = 0.05 [Ω] and E = 0.35 [V ].
We also set the parameter of the Fourier transformation
circuit as M = 10.
Fig. 8 shows the sine-cosine circuit derived from thecircuit in Fig. 7. In this example, for simplify, we consider the basic harmonic component only in the sinecosine circuit and the Fourier transformation circuit. If
we make the circuit for harmonic component by the sameway, we can perform multi-tone analysis.
The simulated results obtained by the conventional
method and the proposed algorithm are shown in Figs. 9
and 10, respectively. The horizontal axis is ω and thevertical axis is the voltage through the nonlinear capacitor C1 . We can see that the step size around theFig. 10. Frequency response obtained by the proposed algorithm (circuit in Fig. 7). Calculationtime: 4.30 [sec].

peaks becomes smaller in the proposed algorithm. We
should mention that the frequency characteristics curve
can be obtained by a single run of the transient analysis of Spice. When we name the peaks of the curve as
peak1, peak2, peak3 and peak4 from the left, the frequencies and the voltages, which should be the exact
values, are summarized in Table 1.
5.2 Example 2
We consider a circuit including
a transistor as shown in Fig. 11. We use the Gummel546

IEEJ Trans. EIS, Vol.131, No.3, 2011

Spice-Oriented Algorithm for Peak Search and Stability Assessment
Table 1. Frequencies and voltages at peaks (circuit in Fig. 7).
ω[rad/s]v1 [mV ]peak1
1.1269523.1

peak2 peak3
3.516 5.2664
862.4 796.3

peak4
6.1803337.3

Fig. 13.

Fig. 11.

Example circuit with strong nonlinearity.

Ampliﬁer circuit including transistor.

Fig. 14. Frequency response obtained by the proposed algorithm (circuit in Fig. 13). Calculationtime: 2.58 [sec].
Table 2. Frequencies and currents at peaks (circuit in Fig. 13).
Fig. 12. Frequency response obtained by the proposed algorithm (circuit in Fig. 11). Calculationtime: 3.69 [sec].

ω[rad/s]
i1 [A]

Poon model for the transistor, whose characteristics are
given as follows;
IS
QB

exp

where
vBC
vBE
1
1−
−
· · · · · · · · · · · · · · · · · · · · · (26)
QB
vAF
vAF

peak3
6.420
3.4450Fig. 14 shows the simulation result of the frequency
response obtained by the proposed algorithm. Although
we consider the basic harmonic component only again,
the curve is quite similar to that reported in (10) . Similar
to the case of the previous example, our algorithm is very
eﬀective to search exact peak values of the frequency response curve. The obtained peaks are summarized inTable 2.
Next, we assess the stability of the solutions of thisstrong nonlinear circuit.

and IS = 2.64 [pA], τ = 100 [ps], Vt = 0.026 [V ], ΦB =0.75, BF = 1940, Cj0 = 24.1 [pF ], VAF = 150 [V ], FC =0.44, BR = 1.
The circuit parameters in Fig. 12 are set as Em = 0.05,
V1 = 10 [V ], V2 = −10 [V ], R1 = 1.5 [kΩ], R2 = 2 [kΩ],
L = 20 [μH], C1 = 50 [pH] and C2 = 1 [μH].
The simulated results obtained by the proposed algorithm is shown in Fig. 11. We can see that the step
電学論 C，131 巻 3 号，2011 年

peak2
6.151
2.5000size around the peaks becomes smaller as the previous example. The peak value Vout = 4.164 [V ] forω = 2.54 [rad/sec] is obtained.
5.3 Example 3
Fig. 13 shows the third example
circuit including strong nonlinear capacitors. This circuit is the same one which has been investigated in (10) .
In order to conﬁrm the correctness of our algorithm, we
use the same circuit parameters as (10) ; Em = 0.35 [V ],α1 = 2.4, β1 = 12.0, α2 = 1.0, β2 = 5.0, R1 = 0.1 [Ω],R2 = R3 = 0.02 [Ω], L1 = 0.2 [H], L2 = 0.8 [H],
L3 = 0.10101 [H],⎧⎨ e = Em sin ωt
C1 ; v1 = α1 q1 + β1 q1 3 · · · · · · · · · · · · · · · · (27)⎩
C2 ; v2 = α2 q2 + β2 q2 3

vBE
−1
Vt
1
vBE
1
−IS
+
exp
−1
QB
BR
Vt
· · · · · · · · · · · · · · · · · · · · ·(24)
vBE
IS
iB =
exp
−1
BF
Vt
vBE
IS
−
exp
− 1 · · · · · · · · · · · · · ·(25)
BR
Vt
iC =

peak1
5.067
2.4545547

Table 3.
λ1
λ2
λ3
λ4
λ5

Fig. 15.

ω = 4.570−0.622 + j0.554−0.622 − j0.554
0.495 + j0.625
0.495 − j0.625
0.8394

Eigenvalues of Φ(T ).
ω = 4.927
0.726 + j0.825
0.726 − j0.825
0.373 + j0.499
0.373 − j0.499
0.8501

ω = 5.560
2.530
0.545 + j0.638
0.545 − j0.638
0.86620.2862the variational circuit. Note that we have to carry out
the transient analysis ﬁve times with diﬀerent initial conditions to obtain all the components.
After obtaining the full matrix Φ(T ), we calculate theeigenvalues by using MATLAB and assess the stabilityof the periodic solutions.
We assess the stability of the solutions at (ω, i1 ) =(4.570, 1.902), (4.927, 0.628) and (5.560, 2.000) as examples. These points are marked in Fig. 14.
First, we show the elementary matrix for the case ofω = 4.570 [rad/s] (Calculation time: 0.72 [sec]).

Variational circuit of Fig. 13.

The circuit equation can be written as⎧
di1⎪
Em sin ωt = R1 i1 + L1
+ αq1 + βq1 3
⎪
⎪
⎪
dt
⎪
⎪
⎪
⎪
di2
⎪
⎪
+ αq2 + βq2 3⎪ αq1 + βq1 3 = R2 i2 + L2
⎪
dt
⎪
⎪
⎨
di3
αq2 + βq2 3 = R3 i3 + L3
⎪
dt
⎪
⎪
⎪
⎪
dq
1
⎪
⎪
= i1 − i2
⎪
⎪
⎪
dt
⎪
⎪
⎪
⎪
⎩ dq2 = i − i .
23dt· · · · · · · · · · · · · · · · · · · · · (28)0.475 0.222
 0.732 0.515
Φ(T ) =  −0.366 0.096
 −2.755 0.630−1.246 −0.551−1.124
1.219
0.813
1.504
6.490

0.130
−0.131−0.021−0.4960.3054

0.066
0.075 

−0.135 
0.345 −0.723

Second, for the case of ω = 4.927 [rad/s] (Calculationtime: 0.95 [sec]).
0.480 0.104
 0.407 0.827
Φ(T ) =  0.012 −0.081
 1.833 −0.026
0.485 −0.407

If we write the variables as periodic solutions with small
variations:
ik = ik0 + Δik· · · · · · · · · · · · · · · · · · · · · · · · · (29)qk = qk0 + Δqk

0.041−0.6401.360−3.0772.228−0.178
0.043
0.121
0.392
−0.045

−0.0160.301 
−0.273 
0.317 −0.010Lastly, for the case of ω = 5.560 [rad/s] (Calculationtime: 0.53 [sec]).

we obtain the following variational equations as
⎧
dΔi1
⎪
⎪
0 = R1 Δi1 + L1
+ αΔq1 + β3q10 2 Δq1
⎪
⎪
dt
⎪
⎪
⎪
⎪ αΔq + β3q 2 Δq
⎪
1
10
1
⎪
⎪
⎪
⎪
⎪
⎪
dΔi2
⎪
⎪
+ αΔq2 + β3q20 2 Δq2= R2 Δi2 + L2
⎪
⎨
dt
dΔi3
⎪
⎪
αΔq2 + β3q20 2 Δq2 = R3 Δi3 + L3
⎪
⎪
dt
⎪
⎪
⎪
⎪
dΔq
⎪
1
⎪
⎪
= Δi1 − Δi2
⎪
⎪
dt
⎪
⎪
⎪
⎪
⎪
⎩ dΔq2 = Δi2 − Δi3dt· · · · · · · · · · · · · · · · · · · · · (30)1.450 −0.083
−0.135 1.001
Φ(T ) =  −0.486 −0.037
 −2.186 0.990
1.537 −0.312−0.611−0.6302.084−2.772−0.7070.154−0.2040.035−0.0720.124



−0.0580.289 

−0.219 
0.4769 
0.309

Table 3 shows the calculated eigenvalues of Φ(T ) for
the three cases.
We can see that the point (ω, i1 ) = (4.570, 1.902) isstable, because all of eigenvalues satisfy |λ| < 1. However, for the other two points, (ω, i1 ) = (4.927, 0.628)
and (5.560, 2.000), the solutions are unstable, becausesome eigenvalues do not satisfy |λ| < 1. These results
agree with the results obtained in (10) . Namely, we cansay that our algorithm gives the correct results withmuch simpler Spice-oriented algorithm.

where we neglect higher-order small terms. From theseequations, we can make the variational circuit of Fig. 13
as shown in Fig. 15. In Fig. 15, q10 and q20 are the steadyperiodic solution and are calculated as

6. Discussions
In 5.3, we assess the stability with all of variables.
However, the states of the circuit elements should bein interdependency. From this reason, in this section,
we try to assess the stability of the solutions with only
limited variables in this section.
For the circuit in Fig. 13, we choose only three variables out of ﬁve state variables and perform the algorithm; namely run transient analysis of Spice for thevariational circuits only three times with three sets of

q0 = Qc cos ωt + Qs sin ωt · · · · · · · · · · · · · · · · · · (31)where Qc and Qs are given by the sine-cosine circuitobtained from the circuit in Fig. 13.
We perform the transient analysis of this circuit justfor one period and obtain the components of the elementary matrix Φ(T ) as the values of the state variables of548

IEEJ Trans. EIS, Vol.131, No.3, 2011

Spice-Oriented Algorithm for Peak Search and Stability Assessment
Eigenvalues when we choose (i2 , i3 , q1 ).

Table 4.
λ1
λ2
λ3

Table 5.
λ1
λ2
λ3

ω = 4.570
0.975
0.260
-0.402

ω = 4.927
1.09260.743 + j0.358
0.743 − j0.358

Acknowledgment
The authors are grateful to Ms. Akiko Kusaka forher work on this topic during 2008-2009. This work waspartly supported by CASIO Science Promotion Foundation.
(Manuscript received July 29, 2010,
revised Nov. 6, 2010)

ω = 5.560
0.212
0.756
2.0449Eigenvalues when we choose (i3 , q1 , q2 ).
ω = 4.570
0.187 + j0.652
0.187 − j0.652-0.780ω = 4.927
0.738 + j0.759
0.738 − j0.759
0.265

ω = 5.560
2.141
0.286
-0.106

References(1)

initial values. The simpliﬁed elementary matrix becomes3 × 3 and we assess the stabilities of the solutions byonly three eigenvalues.
Firstly, for the solution at (ω, i1 ) = (4.570, 1.902),
which is stable, we tried to assess the solution using different ten patterns of variables. Eight patterns out of ten
gave the correct result (stable), while two gave the wrong
result (unstable). Secondly, for (ω, i1 ) = (4.927, 0.628),
which is unstable, only four patterns gave the correct
result. Lastly, for (ω, i1 ) = (5.560, 2.000), which is also
unstable, nine patterns gave the correct result.
Although we do not have obtained any criteria, correct
results are always obtained when we choose (i2 , i3 , q1 )or (i3 , q1 , q2 ). Table 4 and Table 5 show the obtained
eigenvalues when we choose (i2 , i3 , q1 ) and (i3 , q1 , q2 ),respectively.
If we analyze a circuit with weak nonlinearity and withrelatively uniform parameters, we may assess the stability with limited variables and we may give some guidelines to choose the necessary variables. However, if acircuit has strong nonlinearity or nonuniform parameters, it must be very diﬃcult to do that. For the next
step of our research, we need to work on the necessary
number and important nodes of the variables to be chosen under a certain conditions on the circuit.

(2)
(3)
(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

7. Conclusions
We have proposed a Spice-oriented peak search andstability assessment algorithm for nonlinear circuits. By
combining the sine-cosine circuit, the Fourier transformation circuit, the solution tracing circuit with the differentiator and the nonlinear limiter, we can search theexact peaks eﬃciently for nonlinear circuits. By introducing Floquet theory, we can perform the stability assessment of the solutions. The simulation results forexample circuits showed that the proposed algorithm
worked well. We also considered the case that we chose
only limited numbers of the valuables for the stability
assessment.
The proposed method cannot be applied to oscillator
analysis directly, because we have to input the oscillation frequency to the circuit. However, if we combine
the proposed method with the algorithm proposed in ourprevious research (4), which can obtain oscillation frequencies using Newton homotopy circuit, the proposed
method would be extended to the oscillator analysis.
As a future research for our proposed method, we would
like to assess stability for printed circuit boards using
our proposed method that includes limited conditions.
電学論 C，131 巻 3 号，2011 年

P. Wambacq and W. Sansen: Distortion Analysis of Analog
Integrated Circuits, Kluwer Academic Publishers (1979)
D. O. Pederson and K. Mayaram: Analog Integrated Circuits
for Communication, Kluwer Academic Publishers (1991)
A. H. Nayfeh and D. T. Mook: Nonlinear Oscillations, WileyInter Science (1998)
Y. Yamagami, Y. Nishio, and A. Ushida: “Analysis of Reactance Oscillators Having Multi-Mode Oscillations”, IEICE
Trans. Fundamentals, Vol.E89-A, No.3 (Mar. 2006)
J. Kawata, Y. Taniguchi, M. Oda, Y. Yamagami, Y. Nishio,
and A. Ushida: “SPICE-Oriented Frequency-Domain Analysis
of Nonlinear Electronic Circuits”, IEICE Trans. Fundamentals, Vol.E90-A, No.2 (Feb. 2007)
T. Kinouchi, Y. Yamagami, Y. Nishio, and A. Ushida: “SpiceOriented Harmonic Balance Volterra Series Methods”, Proc.
of NOLTA’07, pp.513–516 (Sep. 2007)
F. L. Traversa, F. Bonani, and S. D. Guerrieri: “A FrequencyDomain Approach to the Analysis of Stability and Bifurcations in Nonlinear Systems Described by DiﬀerentialAlgebraic Equations”, Int. J. Circuit Theory and Appl.,
Vol.36, pp.421–439 (2008)
Y. Inoue: “DC Analysis of Nonlinear Circuits Using SolutionTracing Circuits”, Trans. IEICE, Vol.J74-A, pp.1647–1655(1991)
A. Ushida, Y. Yamagami, I. Kinouchi, Y. Nishio, and Y. Inoue: “An Eﬃcient Algorithm for Finding Multiple DC Solutions Based on the SPICE-Oriented Newton Homotopy
Method”, IEEE Trans. Computer-Aided Design of Integrated
Circuits Syst., Vol.21, No.3, pp.337–348 (Mar. 2002)
E. Ikeno and A. Ushida: “The Arc-Length Method for theComputation of Characteristic Curves”, IEEE Trans. Circuits
Syst., Vol.23, pp.181–183 (1976)
A. Ushida and L. O. Chua: “Tracing Solution Curves of Nonlinear Equations with Sharp Turning Points”, Int. J. Circuit
Theory and Appl., Vol.12, pp.1–21 (1984)

Hiroshige Kataoka (Student Member) was born in Osaka,
Japan, on November 25, 1987. He received
the B.E. degree from Tokushima University,
Tokushima, Japan, in 2009. His research interests include computer-aided simulations of
nonlinear circuits.

Yoshihiro Yamagami (Non-member) graduated from the
Department of Electrical Engineering, Anan
Technical College, Tokushima, Japan, in 1982.
He joined the Department of Electric Engineering of Tokushima University, Tokushima,
Japan, in 1984. He received the Ph.D. degrees
in Electric Engineering from Tokushima University, Tokushima, Japan, in 1996. His research interests include computer-aided simulations of nonlinear circuits and high-speed
interconnect circuits.

549

Yoshifumi Nishio (Non-member) received B.E., M.E., andPh.D. degrees in electrical engineering from
Keio University, Yokohama Japan, in 1988,
1990, and 1993, respectively. In 1993, hejoined the Department of Electrical and Electronic Engineering at Tokushima University,
Tokushima Japan, where he is currently a Professor. From May 2000 he spent a year in the
Laboratory of Nonlinear Systems (LANOS) at
the Swiss Federal Institute of Technology Lausanne (EPFL) as a Visiting Professor. His research interests include analysis and application of chaos in electrical circuits, analysis of synchronization in coupled oscillatory circuits, developmentof analyzing methods for nonlinear circuits, theory and applicationof cellular neural networks, and neural network architecture.

Akio Ushida (Non-member) received the B.E. and M.E. degrees in electrical engineering from Tokushima
University, Tokushima, Japan, in 1961 and
1966, respectively, and the Ph.D. degree inelectrical engineering from University of Osaka Prefecture, Osaka, Japan, in 1974. He
had been a Professor in the Department ofElectrical Engineering, Tokushima University,
Tokushima, Japan. From 1974 to 1975, hespent one year as a Visiting Scholar at the Department of Electrical Engineering and Computer Sciences at the
University of California, Berkeley. His current research interests
include numerical methods and computer-aided analysis of nonlinear circuits.

550

IEEJ Trans. EIS, Vol.131, No.3, 2011