2004 RISP International Workshop on Nonlinear
Circuit and Signal Processing (NCSP'04)
Hawaii, USA, Mar. 5-7, 2004.

Finding Complex Conjugate Solutions of Nonlinear Algebraic Equations
Seizo Hagino† , Koji Takata† , Yoshifumi Nishio† and Akio Ushida†††

Tokushima University
2-1 Minami-Josanjima, Tokushima, Japan
Phone:+81-88-656-7470, FAX:+81-88-656-7471
Email: hagino@ee.tokushima-u.ac.jp††
Tokushima Bunri University

abstract: Numerical analysis of nonlinear algebraic
equations is one of the most fundamental and important
problems in scientiﬁc and engineering ﬁelds. In this paper, we propose an eﬃcient SPICE-oriented algorithm
for calculating the multiple solutions and bifurcation
points. Firstly, a given nonlinear equation is replaced by
the equivalent resistive circuit model using ABMs (analog behavior models) of SPICE. Therefore, the DC solutions correspond to those of the nonlinear equations.
We apply the Newton homotopy method to solve theequations, where the homotopy paths are traced by aneﬃcient algorithm based on the arc-length method. To
ﬁnd all the solutions, it is important to ﬁnd the limit
points and pitchfork bifurcation points, because the solution curve in the real space will be disappear at thelimit points, and happens to appear the complex conjugate solutions from the point. Therefore, we propose
an algorithm for calculating the solution curves in thecomplex space starts from the limit point. In such way,
we can ﬁnd out all the solutions in most of the cases.

(1)

There have been proposed many eﬃcient algorithms using SPICE simulators [6]-[8] whose algorithms are basedon the homotopy methods. Although they can eﬃciently
ﬁnd the multiple solutions, it is not known whether they
have gotten all the solutions or not. The homotopy paths
may consist of the multiple independent branches andthe closed loops, whose branches may cross at the points
called pitchfork bifurcation points. Unfortunately, thecurve tracing algorithm based on arc-length method [10]may fail to trace the curve, because the rank of the Jacobian matrix of eq.(1) is reduced to equal or less thann − 1 at the points. On the other hand, since the rank isequal to n at the limit point, we can successfully continue
to trace the homotopy path. The bifurcation points belong to the singular points, so that the curve may turnor cross the other curves there. Therefore, it is important to calculate the bifurcation points to get multiple
solutions or to know the properties of nonlinear equation(1).
In this paper, we show the curve tracing algorithm using SPICE in section 2, and SPICE-oriented algorithms
to ﬁnd the limit and pitchfork bifurcation points in section 3. We also show that the limit point in the real
space also corresponds to the limit point of the complex
space. We found an interesting property such that, bytracing the solution curve in the complex space, we canﬁnd another isolated closed loop in the real space. In
this way, we can calculate the multiple solutions or all
the solution of nonlinear algebraic equations. The interesting illustrative examples are shown section 3 and
4.

Some of them [1]-[2] can eﬃciently ﬁnd all the solutionsof piecewise-linear resistive circuits, where the nonlinear
elements are modeled by the piecewise-linear segments.
Although the interval methods [3]-[4] can also ﬁnd all thesolutions, it is rather time-consuming for large scale systems. The homotopy methods [5] are usefully applied to
calculate the multiple solutions of nonlinear equations.

2. Curve tracing algorithm
Nowadays, SPICE is widely used circuit simulationsof integrated circuits. Since nonlinear devices in ICs
are modeled by many kinds of the functions, SPICE has
all kinds of the controlled sources, and nonlinear functions such as exponential functions, multiplications, additions and subtractions and so on. Furthermore, we
can easily get any function using ABMs(analog behav-

1. Introduction
It is very important to calculate the multiple solutionsof nonlinear equations, because there are many kind ofresistive circuits such as ﬂip-ﬂops and Schmitt-trigger
which have the multiple DC solutions. Unfortunately, awell-known Newton-Raphson method cannot any more
be applied to ﬁnd the multiple solutions. Therefore,
there have been published many papers [1]-[5] for ﬁnding the multiple and/or all the solutions of the following
nonlinear equation:
f1 (v1 , v2 , · · · , vn ) = 0 
f2 (v1 , v2 , · · · , vn ) = 0

fn (v1 , v2 , · · · , vn ) = 0

- 427 -

ior models) in such manner that the desired functions
are programmed with Fortran language. An example ofsingle-output function having the multi-inputs is shownin Fig.1, where the function f (v1 , v2 , . . . , vn ) in the blackbox is realized by writing the function with the input
variables {v1 , v2 , . . . , vn }. On the other hand, a set of
the nonlinear equations given by (1) is realized by theequivalent circuit as shown in Fig.2.
output

input

v1
v2
...
vn

Analog Behavior Model+f(v1,v2,...,vn)
-

ρ

ρ-1
VCCS

-1
STC

.

I 0 i , i=1,2,...nnNonlinear(ρ−1)I0i

resistive circuit

i=1,2,...n.

ρ (=v n+1 )Initial guessvi i=1,2,...n.2 +v= Σ vi
i=1

.ρ2

v

-

ρ

R
1[A]
INTEG

(a)

vccs

vccs

(b)

Figure 3: (a)Equivalence circuit of Newton homotopy method.
(b)Solution curve tracing circuit.

Figure 1: Analog Behavior Model

3. Calculations of bifurcation points
Now, consider an algorithm for calculating the bifur+++
cation points of the nonlinear equations. In this section,
v1 1[Ω]
v2 1[Ω]... vn 1[Ω]we show an eﬃcient SPICE-oriented algorithm. For simplicity, we set vn+1 = ρ and rewrite the relation (3) as
follows:
f1(v1,v2,...vn)+v1 f2(v1,v2,...vn)+v2 fn(v1,v2,...vn)+vn
F(v) = 0
F : Rn+1 → Rn(5)
Figure 2: Equivalence circuit of the nonlinear algebraic equation
Namely, if we ﬁrst assume the voltage at ith resistor It is proved that our algorithm [9] can trace the solution
curves only if the rank of the Jacobian matrix of (5) is(1Ω) voltage vi , we can realize the current source.
n on the whole curve. Since the rank at the limit points
Ii = fi (v1 , v2 , . . . , vn ) − vi , i = 1, 2, . . . , n(2) is equal to “n”, and less than or equal to n − 1 at thepitchfork bifurcation points. Therefore, the limit points
Thus, we have fi (v1 , v2 , . . . , vn ) = 0 with the Kirchhoﬀ’s can be traced by our arc-length method. However, it
voltage law. Let us apply the Newton homotopy method is usually not so easy to calculate the Jacobian matrix
to solve the equation (1). Then, we formulate the fol- and the rank. Thus, let us to obtain an algorithm tocalculate the bifurcation points without estimating thelowing modiﬁed relation:
Jacobian matrix.
(3)
F(v, ρ) = f (v) + (ρ − 1)I0 = 0
For the pitchfork point of (1), we consider a set of thefollowing
modiﬁed equations1 :where f (v) : Rn → Rn given by the relation (1), and
I0 = f (v0 ) is the initial guess. Thus, the relation (2) isF1 (v, h) = F(v) + hI = 0
consisted “n” equations having “n + 1” variables, so that
(6)
F2 (v, ∆v) = F(v + ∆v) − F(v) = 0
the solution satisfying the relation (3) is given by thesolution curves in “(n + 1)” dimensional space. In order
to trace the solution curve with the arc-length method, which have 2n+2 variables and 2n equations if we assume∆v as the additional variables. Since the rank of thewe apply the following relation [9]:Jacobian matrix of F(v) is equal or less than n − 1 [10],K 
if we set the variables ∆v1 and ∆v2 small ﬁxed constants
dvi 2=1(4) in F2 (v, ∆v)), then, eq.(6) is the nonlinear equation withds
“2n” variables and “2n” equations. Hence, we also apply
i=1
the Newton homotopy method as follows:
Observe that the solutions satisfying (1) are obtained
F1 (v, h) + (λ − 1)I10 = 0
at ρ = 1 on the solution curve of eq.(3). Now, we show
F2 (v, ∆v) + (λ − 1)I20 = 0
the equivalent circuit satisfying the Newton homotopy
KK(7)  dvi 2   d∆vi 2  dh 2  dλ 2method (3) and (4). The block diagram is shown by
+
+
+
=1 


ds
ds
ds
ds
Fig.3(a). Some of the “K” node voltages in the equation
i=1
i=3(1) are chosen as the variables of the arc-length method(4) [8]. Thus, the circuit realizing the arc-length method Tracing the homotopy path (7), we can ﬁnd the bifuris shown by Fig.3(b). that is called STC(solution curve cation points at λ = 1 and h = 0. Note that two of the
tracing circuit)[6]. In many practical problems [8], we variables ∆v1 and ∆v2 and the initial guess I20 must be
can get all the DC solutions from a suitable initial guess chosen the ﬁxed suﬃcient small values.
1 At the pitchfork points, the second order term of the Taylor
I0 in equation (3).
expansions are assumed negligible values for the small variations- 428 -Example 3.1. Now, let us calculate bifurcation pointsof the following equation:

4

3

2

2
1

0

0

ρ

-2

-2

-6

1
f (v1 , v2 ) = (v13 − 4v1 − v2 )(− v13 − v2 + 1) = 0
4

-3

-8

(8)

ρ

-1

-4

-10

-5

0

5
s[sec]

10

15

-4
-5

20

0

(a) (2.0,0.0)5s[sec]

1015(b) (-2.0,0.0)

3Thus, we have the following modiﬁed relations:1f1 (v1 , v2 , h) = f (v1 , v2 ) + h = 0
f2 (v1 , v2 ) = f (v1 + ∆v1 , v2 ) − f (v1 , v2 ) = 0f3 (v1 , v2 ) = f (v1 , v2 + ∆v2 ) − f (v1 , v2 ) = 0

0

40

ρ
20

-1
v2

-2

(9)

-4
-5
-5

-20

0

5

10
s[sec]

1520

25(c) (-0.25,0.98)∂f (v1 , v2 )∂f (v1 , v2 )= 0 and=0∂v1∂v2

at the pitchfork points. Hence, the ﬁxed small variables∆v1 and ∆v2 can be set as shown in (9). Applying the
Newton homotopy method, we have
f1 (v1 , v2 , h) + (λ − 1)I01 = 0
f2 (v1 , v2 ) + (λ − 1)I02 = 0
f3 (v1 , v2 ) + (λ − 1)I03 = 0
(10)

2
2
2
2

dv1
dv2
dh
dλ


+
+
+
=1 
ds
ds
ds
ds

0

-3

Note that (8) consists of two variables and one equation,and that it satisﬁes both

60

2

-40
-10

-5

0
v1

5

10

(d) Solution curves of eq.(8)Figure 5:
4. Solution curves of algebraic equations in complex space
We consider the solution curves of the following equation:MN

aij v1i v2j = 0

f (v1 , v2 ) =(12)

i=0 j=0

Assume that the solution curves may consist of several
disjoint loops and/or independent branches. The independent branches with a pitchfork point can be found by
where we chose ∆v1 = ∆v2 = 10−5 . Then, we have the above method in section 3. In this section, we will
3 pitchfork bifurcation points (2.0, 0.0), (−2.0, 0.0) and develop an algorithm to ﬁnd the disjoint loops. Equation(12) satisﬁes the following relations:(−0.25, 0.98) at h = 0, λ = 1 as shown in Fig.4.
∂f (v1 , v2 )∂v2
= 0 or equivalently
=0
∂v1
∂v1

15
10

at the limit points. Therefore, the points can be also
found by solving
f1 (v1 , v2 , h) = f (v1 , v2 ) + h = 0(14)
f2 (v1 , v2 ) = f (v1 + ∆v1 , v2 ) − f (v1 , v2 ) = 0

h

5
ρ=1
0
-5

ρ

-10
-60

-50

(13)

-40

-30

-20
s[sec]

-10

0

10

20

Figure 4: Solution of the pitch fork point
The algorithm is much more simple compared to themethod [11] using the Jacobian matrix of (5), especiallyfor the large scale systems.
Now, consider to calculate the branches starting fromthe pitchfork bifurcation points in (5). We set a small
sphere centered at the bifurcation points v̂, and get theintersections between the branches of (5) and the sphere.
Thus, we haveF(v) = 0 ,K(vi − v̂i )2 = γ

(11)

i=1

for a suﬃcient small positive value γ. The relations (11)are also solved by the Newton homotopy method. The
results for the above example are shown by Fig.5(a),(b)and (c), respectively. Starting from these intersections,
we can trace the solution curves of (5) with arc-length
method [9].

by choosing a suﬃciently small ﬁxed constant ∆v1 .
Thus, the limit point (v̂1 , v̂2 ) can be obtained by solving(14). Note that when the value v2 is increased or decreased from the limit point v̂2 , the solutions v1 of (12)will be given by the complex conjugate values or the realvalues. We consider here the solution curves in the complex space. In this case, we set the variables (a + jb, v2 )instead of (v1 , v2 ) in (12). The directions of the branches
are found by the intersections of a small sphere with thebranches of (12) as follows:{f ((a + jb), v2 )} = 0
{f ((a + jb), v2 )} = 0(15)(v̂1 − a)2 + b2 + (v2 − v̂2 )2 = r2
for a small ”r”. Thus, the solution curves in the complex
space can be found by solving the ﬁrst two equations in(15). Remark that a new solution curve satisfying (12)in real space will be also appeared at another limit pointin the complex space. Thus, the other disjoint loop of(15) in real space can be found in the same way.

- 429 -Example 4.1 Now, consider a simple example as fol- curves are calculated. Thus, the solution curves startlows:ing from the these points are obtained by the arc-length
method [9]. The homotopy paths of the nonlinear equaF (v1 , v2 ) = v14 − v13 − 2v12 − 2v1 − v2 + 4 = 0(16) tion sometimes consist of the disjoint closed loops. To
ﬁnd the multiple loops, we need to trace the solution
For ﬁnding the limit point, we calculate the Jacobian
curves in complex space starting from the limit points.
We also showed an interesting result that the limit point∂F (v1 , v2 )= 4v13 − 3v12 − 4v1 − 2 = 0(17) in real space is also the limit point on the solution curves∂v1
in the complex space. Although it is impossible to prove
The real solution of (17) corresponding to the limit point that all the solutions of the nonlinear algebraic equaof (16)is given by (v̂1 = 1.58, v̂2 = −1.81). Now, set tions are found with our method, the all solutions canv1 = a + jb in (16) for getting the solution curve in be found in many practical problems.
complex space. Then, we have
In the future problems, we need to extend the algo4
3
2
rithm
to large dimensional cases, and also extend to
{F (a + jb, v2 )} = a − a − 2a − 2a + 4− (6a2 − 3a − 2)b2 + b4 − v2 = 0(18) the nonlinear equation instead of the nonlinear algebraic
{F (a + jb, v2 )} = 4a3 − 3a2 − 4a − 2 − (4a − 1)b2 = 0
equations.
References
We assume a small sphere centered at the limit point asfollows:[1] K.Yamamura and T.Ohshima, “Finding all solutions of(19)(a − v̂1 )2 + b2 + (v2 − v̂2 )2 = r2
piecewise-linear resistive circuits using linear programming,”
Substituting (v̂1 = 1.58, v̂2 = −1.81) and r = 0.1 into
(19) and solving (18) together with (19), we have theintersection points (a = 1.58, b = ±0.0245, v2 = −1.87).
Starting from the intersections on the sphere, we canget the solution curve in the complex space as shown in
Fig.6.
Real solution curve
Complex solution curve

IEEE Trans. on Circuits and Systems-I, vol.45, pp.434-445,
1998.
[2] K.Yamamura, M.Sato, O.Nakamura and T.Kumakura, “Aneﬃcient algorithm for ﬁnding all DC solutions of piecewiselinear circuits,” IEICE Trans. Fundamentals, vol.E85-A,
PP.2459-2468, 2002.
[3] L.Kolev and V.Mladenov, “An interval method for ﬁnding
all operating points of nonlinear resistive circuits,” Int. Jour.
Circuit Theory and Applications, vol.18, pp.257-267, 1990.
[4] L.Kolev, “An interval method for global nonlinear analysis,”IEEE Trans. on Circuits and Systems-I, vol.47, pp.675-683,
2000.

υ2

[5] W.I.Zangwill and C.B.Garcia, Pathways to Solutions, Fixed
Points and Equilibria, Prentice-Hall, Inc. Englewood Cliﬀs,N.J., 1981.

20
15
L1

10
5

L2

S4

S2

S1

[6] Y.Inoue, “A practical algorithm for dc operating-point analysis of large scale circuits,”IEICE Trans. Fundamentals,
vol.J77-A, pp.388-398, 1994.

L3

0
-5

S3

-10

3

[7] R.C.Melville, Lj.Trajkovic, S.C.Fang and L.T.Watson, “Artiﬁcial parameter homotopy methods for the dc operating
point problem,” IEEE Trans. on CAD, vol.12, pp.861-877,
1993.

2
-1.5

1
-1

0

-0.5

b

0

0.5

a

-1
1

1.5 -2

Figure 6: Solution curves in the complex space
5. Conclusions and remarks
It is very important to ﬁnd all the solutions of nonlinear equations. In this paper, we have shown Newton
homotopy method for ﬁnding the multiple solutions of
nonlinear equations. The solutions are located on thesolution curves satisfying the modiﬁed equation. The
branches may consist of the disjoint and/or closed loops.
We proposed an eﬃcient SPICE-oriented algorithm forcalculating the bifurcation points such as limit and/or
pitchfork points. After then, the intersections of a small
sphere centered at the bifurcation points with solution[8] A.Ushida, Y.Yamagami,Y.Nishio, I.Kinouchi and Y.Inoue,“An eﬃcient algorithm for ﬁnding multiple DC solutions
based on SPICE-oriented Newton homotopy method,”IEEETrans. on CAD, vol.21, pp.337-348, 2002.
[9] A.Ushida and L.O.Chua, “Tracing solution curves of nonlinear equations with sharp turning points,” Int. Jour. Circuit
Theory and Applications, vol.12, pp.1-21, 1984.
[10] M.Kubicek and M.Marek, Computational Methods in Bifurcation Theory and Dissipative Structure, Springer-Vrlag,
New-York, 1983.
[11] N.Yamamoto, “Bifurcations of solutions of nonlinear equations involving parameters” Theoretical and Applied Mechanics, University of Tokyo press, vol.33, pp.435-444, 1985.

- 430 -